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ABSTRACT 

We examine the detailed physics of the feedback mechanism by relativistic AGN jets interacting 
with a two-phase fractal interstellar medium in the kpc-scale core of galaxies using 29 3D grid-based 
hydrodynamical simulations. The feedback efficiency, as measured by the amount of cloud-dispersal 
generated by the jet-ISM interactions, is sensitive to the maximum size of clouds in the fractal cloud 
distribution but not to their volume filling factor. Feedback ceases to be efficient for Eddington ratios 
Pjct/iodd ;$ 10~ 4 , although systems with large cloud complexes > 50 pc require jets of Eddington ratio 
in excess of 10 -2 to disperse the clouds appreciably. Based on measurements of the bubble expansion 
rates in our simulations we argue that sub-grid AGN prescriptions resulting in negative feedback in 
cosmological simulations without a multi-phase treatment of the ISM are good approximations if the 
volume filling factor of warm phase material is less than 0.1 and the cloud complexes are smaller than 
~ 25 pc. We find that the acceleration of the dense embedded clouds is provided by the ram pressure 
of the high velocity flow through the porous channels of the warm phase, flow that has fully entrained 
the shocked hot-phase gas it has swept up, and is additionally mass-loaded by ablated cloud material. 
This mechanism transfers 10% to 40% of the jet energy to the cold and warm gas, accelerating it 
within a few 10 to 100 Myr to velocities that match those observed in a range of high and low redshift 
radio galaxies hosting powerful radio jets. 

Subject headings: galaxies: evolution - galaxies: formation - galaxies: jets - hydrodynamics - ISM: 
jets and outflows - methods: numerical 



1. INTRODUCTION 

The formation of galaxies is a non-linear, but to some 
degree self-regulatory process; the star-formation effi- 
ciencies of galaxies and the growth rate of the central su- 
permassive black-holes (SMBH) are thought to be mod- 
ified by feedback processes from active galactic nuclei 
(AGN) resulting in a tight correlation between SMBH 
mass an d the bulge stellar velocity dispersion (the M-a 
relation, |FerraresefeMerritt|2000[|Gebhardt et al.|2000| 



Tremaine et al 



2002) 



It is unclear, however, which types of AGN activit y are 
relevant i n regu lating bulge and SMBH growth. The Silk. 
& Rees ([1998) model invokes an energy-driven quasar 
id of low Edd 



wind i 



dingt on rat io, while the models b y |Pabian| 



(1999), |King| ( |2003| ), and |Murray et al.| ( 2005[ ) consider 
opacity-regulated momentum-driven outflows requiring 
Eddington ratios of a few percent. Another possibility 
is that the radiation field in the bulge of galaxies con- 
trols the accretion rate of matter into the centr al regions 
( |Umem ura 2001 Kawakatu &: Umemura|[2002 ). Cosmo- 
logical SPfi and semi-analytic models routinely include 
feedback by powerful radio jets or quasar winds, albeit, of 
necessity, using highly simplified models for the feedback. 
Observationally, ionization diagnostics may not conclu- 
sively distinguish the contributions of radiatively driven 
feedback and feedback driven by je t-ISM interactions 
(Holt et al. 2009 Hayashi et al. 1120121), although i n some 



cases jet-ISM interactions are strongly favoured (Dopita 
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et al. 1997 Nesvadba et al. 20101. Several studies find 
statistical correlations between AGN activit y, outflows 



and the suppression of star-fo rmation (e.g. Schawinski 



et al. 2007 Farrah et al. |2012[ ), but the connection be 



tween^.GN jets and star-formation rem ains ambiguous 



picken et al.||2012| [ Hayashi et al.||2012[ ). 

In cosmolog ical~SPH simulations (e.g. [ Okamoto et al 
2008| |Di Matteo et al.||2008| |Schaye et al. 2010) , grid- 



based simulations {e.g. |Springel] |2Ull[ [ Dubois et al.| 
2012|), and semi-analyt ic models (e.g. |(Jrot.on et al.|200b[ 
Fanidakis et al.||2012|) AGN feedback" 



is found to be a 
necessary ingredient in order to reproduce the observed 
galaxy luminosity function and its evolution with red- 
shift, but the relevant range of powers varies between 
models. Cosmological SPH simulations require energy 
injection rates described by Eddington ratios 77 > 10~ 2 
while some semi-analytic models find that low-powered 
injection of energy with Eddington ratios of 7/ > 10~ 5 is 
sufficient. In both methods there exist a variety of "sub- 
grid" prescriptions to deposit energy yielding different re- 
sults. Neither method resolves or treats the galaxy-scale 
physics of the interaction of the outflows and interstellar 
medium (ISM) adequately, and one of our aims is to pro- 
vide a a robust description of sub-grid feedback physics 
that can be used in future semi-analytic and cosmological 
models. 

Feedback involving mechanical energy input by an 
AGN jet, often termed "radio-mode" feedback, has been 
identified a s a key mechanism to heat the IGM o f the 
cluster (e.g. |Binney fc Tabor|1995j|Soker et al |2001[ ) and 
prevent a runaway build-up of galaxy mass thro ugh fur- 
ther accretion of cooling gas (see|Best et al. 20061 and ref- 



erences therein). Well-studied nearby examples include 
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the H ydra A (|Wise et al.|2007| ), Pers eus A (Fabia n et al. 
2006| ), and M87 in the Virgo cluster ( |Million~eTal.|2U10[ ), 



and the phenomenon is well reproduced in clustcr-sca 



(2011) 



grid-b a sed hydrodynamic simulatio ns by Gaspari et al. 
20121 , IDubois et aT] piToj [20TT| ) and |Teyssier et aT 



Galaxy-scale jet-regulated star-formation ("positive' 
feedback) may be very relevant at higher redsh ifts in gas 
rich galaxies and proto-galactic environments flDe Young) 
T989[ |Bicknell et al.|[2000l |Reuland et al.||2003[ |Rlamer 



et al. 2004; Milcy et al. 2006; Villar-Martm 2007; Vene- 
mans et al.|2007||Miley fc De Breuck|2008HHayasni et al. 



2012p . The case for jet-induced star-formation in disc 



galax i es was made in numerical work as early as [Wood- 
ward| ( 1976 ) and in more recent simu lations by [Fragile 
eTaT^ ( p07jg| and [Gaibler et al | ( pOlT] ). 

In some nearby and high-redshift radio galaxies 
(HzRG) neutral and line-emitting gas is see n outflowing 
at several lOOkms -1 to several lOOOkms -1 ( Gelderman 



fc Whittle 1994; Tadhu nteret al.|2001[ [O^Uea et al.|!M)ij| 
Emonts et al.|2005||Holt et al |2008i| 2011||Morganti et all 

20101 



2005H2007H20T 0; Ncsva cIbTet al.|2006[|20071|20081. 
Lehnert et al. 120111 |Dasyra fc C ombes 2012; Guillard 



et al.|2012||Torresi et al.|2012|). The alignment of the jet 



with outflowing gas ( Pentericci et al.||200l] |Privon et al. 
2008 1, and matching energetics ( |lNesvadba et al. 2006 
2007) suggest that the outflows are driven by the trans 



ter oi energy and momentum from the jet to to the dense 
ISM. This hypothesis is supported by our previous 3D 
hydrodynamical simulatio ns of AGN-jet driven outflows 
( [Wagner fc Bicknell||2011[ WB11 henceforth). 

The question ot how a collimated jet may impart en- 
ergy and momentum isotropically, e.g., to affect the en- 
tire volume defining the bulge of a galaxy, is fre quently 
mentioned ( |De Young||2010| |Ostriker et al.||2010[ ). A re- 
lated problem is the momentum budget associated with 
the dispersion or expulsion of clouds in a galaxy. An im- 
portant feature of AGN jets is that the jet is extremely 
light and that jet and cocoon are highly overpressured 
(under expanded) with respect t o the ambient env iron- 
ment d Begelman fc Cioffi||1989|). Simu l ations by Sax- 



309), 



ton et aITP005| |, [Sutherlan d fc Bicknell| ( [2007] , |Gaibler 



and WB11 ot AGJN jets show how the light, 



overpressured jet inflates a cocoon that drives a quasi- 
spherical energy bubble into the ISM. These simulations 
also showed that isotropization of the injected energy 
is even more effective when the jet encounters inhomo- 
geneities because, by virtue of its lightness - the jet par- 
ticle density is typically 6 to 8 orders of magnitude lower 
than that of ISM clouds - the jet is strongly deflected by 
the inhomogeneities. Additional effects, e.g., jet insta- 
bilities and jet precession increase the isotropy of energy 
deposition, but are not essential. 

In previous work (WB11) we used grid-based hydro- 
dynamic simulations to model jet-ISM interactions and 
quantified the feedback efficiency provided by relativistic 
AGN jets in the core of young, gas-rich radio galaxies. 
The simulated galaxies typically represent either Com- 
pact Steep Spectrum (CSS) or Gigahertz Peaked Spec- 
trum (GPS) sources (Bicknell et al. 1997), which in our 
view are a class of objects experiencing an early phase of 
powerful jet-mediated feedback. In these objects, radio 
source expansion is impeded by the dense multi-phase en- 



vironment of the galaxy core in the early phase of their 
evolution. We concluded that AGN jet feedback in these 
systems is effective in all galaxies for jets with powers 
10 43 - 10 46 ergs -1 if the ratio of jet power to Eddington 
luminosity r\ > 10" 4 . 

A unique feature of these simulations is the treatment 
of the galaxy ISM with a two-point fractal, single point 
log-normal warm-phase distribution (clouds) embedded 
in a hot atmosphere. We determined feedback efficien- 
cies as a function of some of the parameters describing 
this distribution, e.g. the density and filling factor of 
the warm gas. The ISM properties in HzRG are uncer- 
tain; while large reservoirs of molecular gas and HI are 
known to exist, the volume filling factors of the cold and 
warm gas and the typical sizes of clouds are not known. 
WB11 restricted themselves to volume filling factors of 
the clouds of 0.42 and 0.13, which are probably at the 
higher end of the range of "typical" values. Further- 
more, we did not investigate the dependence on maxi- 
mum cloud sizes. 

With the 15 new simulations presented in this paper, 
we have now substantially extended this parameter space 
study to lower filling factors and a variety of maximum 
cloud sizes in the fractal distribution. We also examine 
the acceleration mechanism in more detail, providing an 
explanation for the high mechanical advantage observed 
by WB11. We describe our methods of computation and 
parameter space next in [[2] and Sj3] and compare the rel- 
evant timescales in the problem in ^4[ We present our 
main results in detail in S|5[ In f[6j we compare our sim- 
ulation results with data from a sample of radio galaxies 
with observed outflows. We also discuss other feedback 
criteria and the review the difficulties in modelling cloud 
ablation. We conclude with a summary of the paper in 

m 

2. EQUATIONS AND CODE 

The system of equations describing the relativistic jet 
plasma, hot atmosphe re, and warm clouds in the one 
fluid approximation is ( Landau fc Lifshitz]|1987[ ): 

dD dDu 1 

1 

dt 



8F l 

~dt 

dE 

~dt 



dF l ui 

dF i c 2 
dx % 



dx % 
dp 
dx' 



0; 
0; 



D 



puTV/c 2 
p 2 A(T); E = pwT 2 -p. 



F' 



(1) 



The conserved quantities, D, F J , and E are the labora- 
tory frame fluid density, components of the momentum 
density, total energy density (including the rest mass en- 
ergy density). The variables p, T, A, and u l are pres- 
sure, temperature, cooling rate, and the components of 
the three velocity, respectively. The bulk Lorentz factor 

1/2 

is r = (l — Mju'/c 2 ) . The proper rest frame density 

is p and w — c 2 +P"f/p('~f — 1) the proper rest frame spe- 
cific enthalpy for an ideal polytropic equation of state, 
with index 7. 

We integrate these equations using the publicly avail- 
able, open-source E ulerian Godunov-type code FLASH 
(Fryxell et al. 20001 version 3.2 and its re lativistic hy- 
drodynamics module (|Mignone et al. 2005 ) to which we 



have added code to incorporate radiative cooling and 
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code to advance advected scalars in the relativistic hy- 
drodynamic solver. 

We exploit the adaptive mesh capabilities of FLASH, 
utilizing up to seven levels of grid refinement in a cubi- 
cal simulation domain of lkpc 3 in physical dimensions, 
consisting of 1024 3 cells at a maximum spatial resolution 
1 pc. This is twice the resolution of the simulations by 
WB11 and is necessary in order to capture the fractal 
outlines of clouds for small filling factors and cloud sizes. 
The jet inlet and initial jet width is 20 pc and is resolved 
with at least 10 cells. The formation of clean diamond 
shocks indicates that the jet stream is sufficiently well 
resolved. Note that a restricted o ne parameter scaling of 
physi cal dimensions is possible (Sutherland & Bicknell 



2007) 



Tracer variables distinguish jet material and warm 
phase gas from each other and from the hot phase 
background. We include non-equ ilibrium, optically thin 
atom ic cooling for T > 10 4 K ( Sutherland fc Dopita 
1993 1 and updated solar abundances ( |Asplund et al. 



2005), for which the mean mass per particle, /i 
0.5165. Thermal conduction, photo-evaporation, self- 
gravity, and magnetic fields are not included. 

We do not include a static gravitational potential for 
the SMBH or the bulge because our simulations span a 
spatial range from 1 pc - 1 kpc, within which neither the 
gravitational force due to the SMBH, nor that due to 
the bulge are dynamically important over the timescales 
considered here. For typical SMBH and bulge masses in 
evolving massive galaxies, the SMBH sphere of influence 
extends to a radial distance of order 10 pc, covering only 
a few tens of cells in our simulation domain around the 
base of the jet. This volume is quickly filled with light jet 
plasma, which is practically unaffected by gravity. On 
the kpc scale, the density and pressure profiles of the 
hydrostatic environment in a massive gas-rich spheroidal 
proto-galaxy are fairl y flat under the grav itational influ- 
ence of the bulge (e.g. Capelo et al.|21)l():. and we adopt 
a uniform hot phase distribution characterized by a tem- 
perature of Th = 10 7 K and a density of cither = 0.1 
or rih = 1.0. The gravitational force due to the bulge 
may be neglected from timescale arguments, described 
hi PI 

We ran our simulations on the National Computational 
Infrastructure National Facility (NCI NF) Oracle/Sun 
Constellation Cluster, a high-density integrated system 
of 1492 nodes of Sun X6275 blades, each containing two 
quad-core 2.93 GHz Intel Nehalem cpus, and four inde- 
pendent SUN DS648 Infiniband switches^] We typically 
used 256 to 1024 cpus with 3GB of memory per core to 
complete one simulation within two weeks. 

3. MODEL PARAMETERS, INITIAL CONDITIONS, AND 
BOUNDARY CONDITIONS 

A crucial ingredient in the simulations described in 

fis the two-phase ISM, which consists of a warm 
- 10 4 K) phase and a hot (T - 10 7 K) phase. In 
particular, we are concerned with the effect of the jet 
plasma on the state and dynamics of the warm phase 
material. We have, therefore, extended our studies of 
parameters related to the warm-phase and identified the 

3 For details of the system spec ifications see 
http://nf.nci.org.au/facilitics/vayu/hard ware.phpl 



correct physical mechanism that leads to the acceleration 
of the clouds. 

The warm phase ISM density is initialized from a cube 
of random numbers that simultaneously satisfies single- 
point log-normal statistics and two-point fractal statis- 
tics. Let P{p) be the log-normal probability density func- 
tion of the random variable p, representing density: 



P{p) 



1 



2irp 



exp 



-(In p — m)* 
2? 



where 



In 



'In 



A* 



1 



P> 



(2) 



(3) 



and \i and a 1 are the mean and variance of the log-normal 
distribution. 

Let F(k) be the Fourier transform of the spatial den- 
sity distribution, p(r), with k and r as wave vector and 
position vector, respectively. The two-point fractal prop- 
erty is characterized in Fourier space by a power spec- 
trum, D(k), in wave number, k, that obeys a power-law 
with index —5/3 for a Kolmogorov-type spectrum 



D{k) 



J k 2 F{k)F*{k)dfl oc fe" 



5/3 



(4) 



where the integral of the spectral density, F(k)F*(k), is 
over all solid angle, il. 

A cube of random numbers that simultaneously satis- 
fies Eqn. (pi) and Eqn. (fH is gene rated by the method 
outlined in ]Lewis fc Ausun (2002). First, a cube with 
cell values from a Gaussian distribution with mean m 
and standard deviation s is Fourier-transformed and 
apodized by a Kolmogorov power law spectrum in 
wavenumber with index —5/3 and minimum sampling 
wavenumber /c m i n . The minimum sampling wavenum- 
ber is, effectively, the average number of clouds per di- 
mension divided by 2, and it determines the scale of the 
largest fractal structures in the cube relative to the size 
of the cube. For example, if fc m i n = 20 for a cube mapped 
to a domain of extent lkpc, then fc m ; n = 20 kpc -1 and 
the largest structures (clouds) would have extents of 
-Rc.max = l/(2fc m in) = 25 pc. The cube is then trans- 
formed back into real space and exponentiated. Be- 
cause the last step alters the power-law structure in 
Fourier space, the cube is iteratively transformed be- 
tween Fourier space and real space until successive cor- 
rections produce a power-law convergence within 1%. 

To place the fractal cube into the simulation domain 
it is apodized (in real space) by a spherically symmetric 
mean density profile which in the simulations presented 
here is flat with mean warm phase density (n w ). The 
porosity of the warm phase arises by imposing an upper 
temperature cutoff for the existence of clouds at T cr i t = 
3 x 10 4 K, beyond which clouds are deemed thermally 
unstable. No lower temperature limit is enforced, and 
temperatures in the core of clouds may initially be less 
than 100 K. The upper temperature cutoff corresponds 
directly to a lower density cutoff, p crit = /i m p/(fcT cr j t ), if 
the pressure, p, is defined. Here, p m is the mean mass 
per particle of the hot phase. In our simulations the 
clouds are in pressure equilibrium with the surrounding 
hot phase, thus p clit = ^ m n h T h /T a - lt , where n h and T h 
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are the hot phase number density and temperature, re- 
spectively. The filling factor of the warm phase, within 
the hemispherical region of radius 0.5 kpc, in which it is 
distributed, is: 



fv 



P(p)d P 



In 



{(Pcrit/M)V^/M 2 + l}' 



(5) 



The original fractal cube was constructed with /i = 1 and 
a = 5, and after apodization with a spatially uniform 
mean warm phase density distribution, the single point 
density distribution remains lognormal, but with a mean 
[i = (n w ). For an isothermal hot phase distribution, 
whose temperature in this work is fixed at Th = 10 7 , 
Pent /A 1 = (T h /T clit )(n h / (n w )) is constant everywhere. 
The filling factor is, therefore, directly defined by the 
ratio of hot phase density and mean warm phase density, 
W (n w ). 

The clouds embedded in the hot phase remain static 
unless impacted by the jet or jet-blown bubble. The gas 
temperatures in most cells containing warm phase ma- 
terial are initially below 10 4 K and are not subject to 
radiative cooling. For a more detailed description of the 
method to generate the fractal cube and a discussion 
of the choice of statistical parameters for the lognormal 
probability distribution and wavenumber powe r law in- 
dex we refer the reader to the manuscript by |Lewis fc| 
Austin ([2002 1 and the relevant sections and appendixes 



Sutherland fc BicknelT| ( |2007| . 

.'he general setup, initial conditions, and boundary 
conditions used here are identical to those of WB11. 
WB11 performed 14 simulations of AGN jets with pow- 
ers in the range 43 < log(Pj et /ergs -1 ) < 46. The 
choice of warm phase filling factors, fv, of 0.42 and 0.13, 
was relatively high and the maximum cloud size fixed at 

-Rc.max ~ 25 pC (fc min = 20kpC~ 1 ). 

In the 15 new simulations presented here, we explore 
new regions of parameter space with filling factors, fv, 
of 0.052 and 0.027, corresponding to average warm phase 
densities of 150 cm -3 and 100 cm -3 , and k m - m of 40 kpc -1 
and 10 kpc -1 , corresponding to maximum cloud sizes of 
max — 12.5 pc and i? c ,max ~ 50.0 pc. The range in 
jet power and other parameters defining the jet plasma 
remain the same. These jets typically have a density 
contrast of 10 -4 with respect to the ambient hot phase 
and 10 -7 with respect to the embedded clouds. The 
pressure contrast between the jet and the hot phase ISM 
is typically 10 2 - 10 3 . AGN jets are extremely light, 
underexpanded (overpressured) jets. 

The complete list of 29 simulations, including those 
from WB11 are given in Table [2j New runs are marked 
with 

4. TIMESCALES 

It is instructive to compare the timescales and associ- 
ated lengthscales present in this problem. The definition 
and values of the relevant timescales are listed in Table[TJ 

An unimpeded jet typically requires of order 100 kyr to 
cross a domain of 1 kpc, while a jet propagating through 
a clumpy ISM is confined anywhere between 20 kyr and 



1 Myr, depending on jet power and average cloud density. 

In comparison to the confinement time, the bulge free- 
fall time is at least two orders of magnitude larger, justi- 
fying our neglect of gravity in the simulations. A closely 
related timescale, the buoyancy timescale for a jet blown 
bubble is also only important on timescales much longer 
than the simulati on time and on spa tial scales much 
larger than 1 kpc ( Bruggen et al. 2002 ) . 



The cooling time in the hot ISM is also longer than 
the simulation time, but the cooling time in the clouds is 
short enough to affect the computational hydrodynamic 
timestep and the cooling length is not resolved in our 
simulations. The sound crossing time inside clouds is 
much longer than that in the inter-cloud medium, and 
also much longer than the jet confinement time, and we 
use this property to slightly underpressure the clouds (by 
~ 2%) to keep the cloud inte rfaces sharp and static. 

|Wagner fc Bicknell (2011) compared the cloud col- 
lapse and cloud ablation timescales and concluded that, 
while clouds engulfed by the jet-blown bubble experi- 
enced an external pressure enhancement that reduced 
the critical Bonnor-Ebert mass sufficiently to formally 
induce collapse, the comparatively short ablation times 
may destroy clouds before stars can form. Cloud ablation 
is facilitated by the Kelvin-Helmholtz instability, which 
grows over timescales comparable to or shorter than the 
ablation timescale. The cloud crushing time is long in 
comparison. These timescales are included here for com- 
pleteness. 

5. RESULTS 

5.1. Velocities of accelerated clouds and feedback 
efficiency 

We conducted 15 new simulations to study new re- 
gions in the space spanned by parameters that describe 
the distribution of warm phase material in our simula- 
tions as described in Sj3] Figure [T] shows density maps 
of three selected new simulations with lower filling factor 
and differing maximum cloud sizes to those of previous 
simulations. To obtain a three dimensional impression of 
the interactions between the jet and the clouds we show 
a volume render of the density of both components from 
one of our simulations in Figure. [2] The jet plasma is 
textured in bluish green and the clouds in purple. The 
forward shock outlining the jet-blown energy bubble is 
seen in a translucent grey. An oval excavation is made 
in the visualization of the clouds in order to show the jet 
plasma flow within. 

The global evolution of the simulations was described 
by Wagner & Bicknell (2011 1. A key feature of the jet- 
ISM interactions is that whatever the initial narrowness 
of the jet, the jet flow is broadened by the interaction 
with the first cloud. The secondary jet streams flood 
through the porous channels of the two-phase ISM and a 
quasi-spherical jet-driven bubble sweeps over the entire 
bulge region. The feedback operates isotropically, with- 
out depending on the initial width or collimation of the 
jet, and the clouds at all position angles in the galactic 
halo are dispersed to high velocities. 

Let Mbhj tn-p, c, <tt, and Cioo be the black hole mass, 
the proton mass, the speed of light, the Thomson elec- 
tron scattering cross section, and the velocity dispersion 
in units of 100 km s -1 , respectively. We also define the 
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Table 1 

Timescales in Jet-ISM interactions 



Timcscalc 



Definition 



Typical values 



Unimpeded jet crossing time( a ' 

Jet confinement time' b ' 
Bulge free-fall time( c ' 
Buoyancy timescale( d ' 
Cooling time in hot ISM 

Cooling time in clouds at critical temperature' 6 ' 
Cooling time in shocked clouds^' 
Cloud sound crossing time's' 
Inter-cloud sound crossing time'' 1 ' 
Cloud Kelvin-Helmholtz growth time W 
Cloud crushing time '■>' 
Cloud collapse time' k ' 
Cloud ablation time' 1 ' 



1 + 



(1+xK 



0.5L/^/2gV/SC 
k B T h /n h A{T h ) 
kBT CI it /n clit A(T clit ) 
k B T s /n a A(T 3 ) 
2R c /c c 
dch/ch 

(Rc/v ch )((n w ) + n c h)/\/ {nw) n ch 
Rc/vs ~ { n w } Rc /n ch v ch 



6 kyr - 300 kyr 

20 kyr - 1 Myr 
120 Myr 
45 Myr 
8.6 Myr 

15.5 yr 

14.6 yr 
6.6 Myr 

66 kyr 

24 kyr 
2.4 Myr 
3.8 Myr 

40 kyr 



' a ' Time required for th e jet head to cross the L = 1 kpc domain if no clouds impeded its progress (see 
e.g. |Safouris et al.|2008| . The variables V and ( denote the lorentz factor and the ratio of jet density to 
ambient gas density, and \ = (7 — l)pj ct c 2 /7Pj ct is the proper density parameter. The lower and upper 



10 



|46 erg s 1 



0.1 cm 3 and R ct 



10 



i43 erg s 1 



values correspond to the cases for which Pj c 
rift = 1.0 cm~ 3 , respectively. 
( b ' Time required for the jet to cross the L = 1 kpc domain in the presence of clouds impeding its 
progress. This is effectively equivalent to the duration of the simulation. 
< c ' Ph = 1.0/w 

( d ' V, S, C, and g are volume and cr oss-section of the b uoyant bubble, the drag coefficient, and the 
;ravitational accelerat ion, respectively j Birzan et al.|2004| . Here we choose V/S = 0.5 kpc and C = 0.75 
Churazov et al.|2001fr 

' n cr ; t is the critical number density corresponding to p cr it • 
( f ' T s and n s are the postshock temperature and particle number density, respectively, of the shock 
propagating into the cloud. A shock speed of 100 km s -1 and preshock conditions of n = n cr it and 
T = T cr it were assumed. 

( g ' c c is the average sound speed in a cloud with average temperature 1000 K. 

( h ' Ch is the sound speed of the hot phase and d c h ~ 2R C is the inter-cloud distance. 

W The growth timescale of the Kelvin-Helmholtz instability at the interface between a cloud and its 

ambient hot flow. jj c j, ~ 10 5 kms" 1 and n c ^ ~ 0.1 cm -3 are the velocity and particle number density of 

the flow through the inter-cloud channels. This is the upper limit corresponding to the lowest excited 

mode. 

0' v s is the speed of the shock propagating into the cloud. 
( k ) This is equivalent to the cloud free-fall time. 

^abl ~ 600 kms -1 is the ablation speed. It corresponds to the channel speed observed in our simula- 
tions. 



ratio of jet power to Eddington luminosity (the Edding- 
ton ratio) to be 77 = Pj e t/i e dd, and 4> Wl p, and v r as 
the warm phase tracer (mass fraction in a cell), density, 
and radial velocity, respectively. A convenient measure 
of the efficiency of feedback is the density-averaged ra- 
dial outflow velocity, (v r>w ) = J2i=i <Pw P v r/J2i=i <l>w P, 
relative to the velocity dispersi on of a galaxy's bulge as 



predi cted by the M-a relation (Silk & Rees 1998 King 
2005 1 . Defining the black hole mass m terms of the Ed 
dington ratio Mbh = ^Gm v P; e tc/riaT and using the 



M-a relation found by Tremaine et al. ( 2002 ) we express 
the velocity dispersion as 



CT ioo 



(G) 



where -Pj e t,45 is the jet power in units of 45 erg s . When 
(v r ,w) > &> the jet-ISM interactions result in sufficient 
feedback of momentum and energy to establish a highly 
dispersed distribution of cold and warm gas within the 
core of the galaxy. The advantage of scaling the jet power 
by the Eddington luminosity is that it allows us to re- 
late a given simulation to the conditions for feedback set 
by the M-a relation in a galaxy with a SMBH mass 
according to the value 77. We may scale the jet power 



arbitrarily because the simulations do not depend on the 
gravitational field for a particular Mbh or bulge mass, 
M*. 

The maximum values of (iv,t«) during the run as a func- 
tion of jet power for all 28 simulations that include the 
warm phase are shown in Fig. [3j which updates Fig. 5 
in WB11. Points of constant hot phase density and fill- 
ing factor are connected along increasing jet power with 
lines of specific colors and style according to the legend. 
The slanted grey lines in Fig. [3] represent the loci of the 
velocity dispersion along lines of constant 77 (as deter- 
mined by Eqn. The dashed grey lines in the figure 
represent a different M-a relation of Mbh oc a 5 , which 
is mo r e in a greement with the recent res ults by|Graham| 
et al. (2011 ), especially for core galaxies ( Graham|2012p . 

The locus in this case is ctioo = 1.2r/ _1 / 5 Pjg{ 5 45 . Using 
either relation, one may then compare the points of the 
simulations with values for the velocity dispersion pre- 
dicted by the M-a relation for a given value of 77. If a 
point lies above an isoline for 77, then feedback by a jet of 
that power Pj ct in a galaxy with Eddington limit Pjet/v 
is effective. Conversely if a point lies below an isoline for 
77, then feedback is not effective. Equivalently, the point 



6 Wagner, Bicknell, & Umemura 



Table 2 

Simulation parameters 



New 


Simulation 


log Pjet (a) 




PlSM/fc (c) 
(cm- 3 K) 








«c,max (g) 








(erg) 


(cm 3 ) 


( cm -3 ) 




(kpc" 1 ) 


(PC) 


(1O 9 M ) 




A 


15 


0.1 




10° 














B 


16 


1.0 


10 7 


1000 


0.42 


20 


25.0 


16 




B' 


46 


1.0 


10' 


300 


0.13 


20 


25.0 


3.2 


► 


B" 


46 


1.0 


,-,7 

10' 


150 


0.052 


20 


25.0 


0.29 


► 


B'" 


46 


1.0 


40 7 


100 


0.027 


20 


25.0 


0.15 




C 


46 


0.1 


40 6 


100 


0.42 


20 


25.0 


1.6 




C 


46 


0.1 


10° 


30 


0.13 


20 


25.0 


0.32 




D 


45 


1.0 


10 7 


1000 


0.42 


20 


25.0 


16 


► 


D/io 


15 


1.0 


10' 


300 


0.13 


10 


50.0 


3.2 




D' 


15 


1.0 


10' 


300 


0.13 


20 


25.0 


3.2 


► 


D i'o 


15 


1.0 


10' 


150 


0.052 


10 


50.0 


0.29 


► 


D" 


15 


1.0 


40 7 


150 


0.052 


20 


25.0 


0.29 


► 


TV' 


15 


1.0 


10 7 


150 


0.052 


10 


12.5 


0.29 


► 


^10 


45 


1.0 


10' 


100 


0.027 


10 


50.0 


0.15 


► 


D'" 


45 


1.0 


40 7 


100 


0.027 


20 


25.0 


0.15 




E 


45 


0.1 


10° 


100 


0.42 


20 


25.0 


1.6 




E' 


45 


0.1 


10° 


30 


0.13 


20 


25.0 


0.32 




F 


44 


0.1 


1Q 6 


100 


0.42 


20 


25.0 


1.6 




F' 


44 


0.1 


40 6 


30 


0.13 


20 


25.0 


0.32 




G 


44 


1.0 


10 7 


1000 


0.42 


20 


25.0 


16 




G' 


44 


1.0 


10 7 


300 


0.13 


20 


25.0 


3.2 


► 




44 


1.0 


10 7 


150 


0.052 


10 


50.0 


0.29 


► 


G" 


44 


1.0 


10 7 


150 


0.052 


20 


25.0 


0.29 


► 


p" 


44 


1.0 


10 7 


150 


0.052 


10 


12.5 


0.29 


► 


G'" 


44 


1.0 


10 7 


100 


0.027 


20 


25.0 


0.15 




H 


43 


0.1 


10° 


100 


0.42 


20 


25.0 


1.6 


► 


H' 


43 


0.1 


10 6 


30 


0.13 


20 


25.0 


1.6 


► 


I" 


43 


1.0 


10 7 


150 


0.052 


20 


25.0 


0.29 


► 


I'" 


43 


1.0 


40 7 


100 


0.027 


20 


25.0 


0.15 



Note. — Runs with run labels containing the same letter are runs with the same jet power, Pj e t, and hot phase density, n^. Runs 
labeled with single, double, or triple primed (" '") letters denote lower filling factor counterparts to runs with less number of primes. All 
runs, other than those whose run label contains the value of A: m i n in the subscript, were performed with fc m i n = 20 kpc . 

( a ) Jet power. 

( b ' Density of hot phase. 

( c ' p/k of both hot and warm phases. 

( d ) Average density of warm phase. 

( e ) Volume filling factor of warm phase. 
W Minimum sampling wave number. 

( g ) Maximum cloud size. 

( h ) Total mass in warm phase. 

velocity of clouds, (v r ^ w ), or equivalently, the critical Ed- 
dington ratio of the jets, ?7 cr it, depends weakly on filling 
factor, but strongly on the maximum size of clouds in 
the galaxy bulge. The overall lower limit f7 cr i t 
for efficient feedback found by WB11 is only slightly 
reduced for galaxies containing small cloud complexes 
(-R c ,max ^ 10 pc, fc min = 40 kpc -1 ) but jets with Edding- 
ton ratios of 77 crrt = 10 -2 - 10 -1 are required if cloud 
complexes are large (i? c ,max ^ 50 pc, /c m i n = 10 kpc -1 ). 

WB11 found that the jet-ISM interactions, despite 
the porosity of clouds and the radiative losses of shock- 
impacted clouds, exhibit a high mechanical advantage, 
meaning that substantial momentum transfer from the 
jet to the clouds occurred through the energy injected 
by the jet. We define the mechanical advantage in our 
simulations at a given time as the ratio of the total ra- 
dial outward momentum carried by clouds to the total 
momentum delivered by the jet up to that time. Fig- 
ure [4] shows the curves for the mechanical advantage as a 
function of time for all 28 simulations including a warm 



itself marks a critical value of 77, 7y cr i t = (J 3 jot/iodd)crit ) 
below which feedback ceases to be effective in galaxy with 
a SMBH of mass M B h = ^Gm p Pj et c/r]a T . 

As observed in previous simulations, the velocities at- 
tained by clou ds match those observed of outflows in 
radio galaxies QMorganti et al.||2005| |Holt et al.||2008l 
Nesvadba et ai.||M8[ MMI |Lehnert et "aTlpnTTj 

Uasyra Combes|201l| |Guillard e t al.|2012||Torresi et al.| 
2012p . The dense cores of the clouds in our simulations 
are accelerated to a few 100 km s -1 , while the diffuse ab- 
lated material is accelerated to several 1000 km s -1 . We 
discovered that the feedback efficiency of the relativistic 
jet on the warm phase ISM increases with increasing jet 
power, decreasing mean ISM density, and increasing fill- 
ing factor, although only two values for the filling factor, 
fv = 0.42, and /y = 0.13, were studied. 

Within the new range of parameter space, the main 
conclusions reached in WB11 remain valid; feedback is 
effective in systems in which the jet power is in the range 
Pjet = 10 43 - 10 46 ergs -1 and ?y > 77 cr it- Furthermore, we 
find that, the maximum density weighted radial outflow 
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Figure 1. Logarithmic density maps (in units of cm -3 ) of selected new simulations. The domain extents in each panel are 1 kpc X 1 kpc. 
The left column of panels show a face on view of initial the warm gas distribution. The center and right columns of panels show midplane 
slices at an advanced stage of the simulations for 2 = (reflected about x = 0) and y = 0, respectively. Top row: Run D'", a very low 



filling factor run (fy = 0.027); Middle row: Run D' 10 , maximum cloud sizes of R c 



sizes of -R c , max = 10 pc. See the electronic edition of the Journal for a color version of this figure. 



50 pc; Bottom row: Run DV Q , maximum cloud 



phase. For all simulations the mechanical advantage is 
much greater than unity. Most curves fall closely on top 
of each other along a narrow band up to at least 1 Myr. 

The high mechanical advantage generally leads to a 
high fraction of jet energy transferred to kinetic energy 
of the warm phase. In Fig. [5] we show the evolution of 
the ratio of kinetic energy in clouds to injected jet en- 
ergy, -Ekin,w/-Pjeti) as a function of t for all 28 simulations 
including a warm phase. In all cases that fraction is high, 
reaching ~ 0.1 - 0.4, with details depending on jet power 
and ISM properties. The details of the dependence on 
feedback efficiency on ISM properties are given in the 
next two sections and the physics of how the high me- 
chanical advantage is sustained and the energy transfer 



occurs are investigated in §5.6| 

5.2. Dependence on filling factor 

Figure [6] shows the maximum values of (v rM ) reached 
in the simulations as a function of fy . The markers de- 
note simulations with equal values of Pj C t/n-h, as indi- 
cated in the legend. The lines of a given color connect 
simulations of equal power, also indicated by the label 
letter, and the line color indicates the hot phase density. 
Apart from the cases of different k nim in the D-series of 
runs, simulations grouped by connected lines, therefore, 
also indicate runs with equal values of P^t/nh- 

In general, the dependence of (v r>w ) on filling factor 
is weak and non-monotonic. In the B-series, D-series, 
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400 pc 



67 kyr 
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296 pc 



125 kyr 




505 pc 



15 B kyr 



187 kyr 
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Figure 2. Volume render of the density of the jet plasma and clouds for run D'. The jet plasma is textured in bluish green and the 
clouds in purple. The forward shock outlining the jet-blown energy bubble is seen in translucent grey. An oval excavation is made in the 
visualization of the clouds in order to show the jet plasma flow within. The view moves outward from the core of the galaxies as the bubble 
of jet plasma expands, and the physical (projected) size is indicated by scale bars on the bottom right in each panel. The simulation data 
is reflected about x = and the left side is rotated by 180° about the jet axis to show a back view of the simulation. See the electronic 
edition of the Journal for a color version and an mpeg animation of this figure. 
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Figure 3. Maximum mean radial velocity of clouds, (v r>m ), 
against jet power for the simulations B — I'" of Table [2] The solid 
and dashed grey lines are loci of constant rj, the ratio of jet power 
to Eddington luminosity, for M— a relations with powers 4 and 5, 
respectively. The line colors indicate different hot-phase densities 
and the line styles represent different filling factors, as indicated 
in the legend. See the electronic edition of the Journal for a color 
version of this figure. 

Dio-series, G-series, and I-series of the simulations (see 
table [2l for nomenclature), we observe that for fy > 0.1, 
lower Tilling factors decrease the feedback efficiency, while 
f° r fv 0.1, lower filling factors increase the feedback 
efficiency. The reason for the weak dependence and the 
non-monotonicity is the competing effects of the cloud 
ablation rate and jet plasma confinement time. On the 
one hand, smaller filling factors increase the volume avail- 
able for the jet plasma to flood through, and thereby re- 
duce the confinement time, which reduces the impulse 
delivered to the clouds over the confinement time. On 
the other hand, smaller filling factors increase the mass 
of ablated material relative to the total mass of a cloud, 
because the ablation rate is proportional to the cloud 
surface area and the mass is proportional to the cloud 
volume, which decreases faster than the former for de- 
creasing filling factor. When lowering the filling factor 
in the range fy > 0.1 the effect of reduced plasma con- 
finement time dominates over the effect of increased frac- 
tional cloud ablation and results in lower mass-averaged 
outflow velocities. In the range fy < 0.1, the increased 
cloud ablation rate dominates over the reduced plasma 
confinement time when reducing fy, leading to higher 
mass-averaged outflow velocities. Over the range of val- 
ues for the filling factor studied here, these two effects 
counteract one another, and the dependence of (v r ^ w ) on 
fv for constant Pjet/nh remains weak. 
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Figure 4. Mechanical advantage versus time for all 28 runs in- 
cluding clouds. The top, middle, and bottom panels shows runs 
for which f v = 0.42, 0.13, and 0.052 or 0.027, respectively. The 
mechanical advantage here is defined as the total outward radial 
momentum of clouds at time t divided by the total momentum 
delivered by the jet up to time t. The mechanical advantage in 
all simulations 2> 1 indicating strong momentum coupling in the 
energy-driven regime. See the electronic edition of the Journal for 
a color version of this figure. 

The mechanical advantage (Fig|4| is slightly reduced 
for systems with lower filling factor down to fy = 0.027, 
but the dependence of the efficiency of transfer of jet 
energy to kinetic energy of the warm phase (Fig[5J) 
on warm-phase filling factor parallels the weak (non- 
monotonic) dependence of the maximum outflow velocity 
on filling factor. 

Note that, by reducing the filling factor, we are also 
reducing the total mass of the warm phase, in contrast 
to this, we may keep the total mass and filling factor the 
same but change the maximum size of clouds by varying 
fcmin- The results for this are shown next. 
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t(kyr) 

Figure 5. Fraction of jet energy going into kinetic energy of the 
warm phase as a function of time for all 28 simulations containing 
a warm phase. The top, middle, and bottom panels shows runs 
for which f v = 0.42, 0.13, and 0.052 or 0.027, respectively. For all 
simulations, 0.4 > S w ,km/-Pj e t > 0.1 although the maxima and the 
time taken to reach the maxima depend on the jet power and ISM 
parameters. See the electronic edition of the Journal for a color 
version of this figure. 

5.3. Dependence on maximum cloud sizes 

Let us look at the D-series of runs, for which we have 
varied the maximum size of clouds by varying k m - ln . The 
values of fc m in are denoted by the subscript of the run 
labels in Fig. [3j 

In Fig. m we plot the sequences D" , D", D4Q, and 
G" , G", TJ40 against fc m i n with black markers, labels, 
and connected by a line. The grey, unlabelled markers 
are other runs with varying filling factors but otherwise 
identical parameters. The sequences in this figure, but 
also those in both Figs. [3] and [6j show that the depen- 
dence of mean velocity on the maximum cloud size in a 
simulation is very strong, and that is much stronger than 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 

fv 

Figure 6. Maximum mean radial velocity of clouds, {v ryW ), ver- 
sus cloud volume filling factor for the simulations B - I'" of Ta- 
ble [2] The line colors indicate different hot-phase densities and the 
marker styles group simulations with equal values of Pj^/n^, as 
indicated in the legend. See the electronic edition of the Journal 
for a color version of this figure. 

the dependence on filling factor. 

By halving the size-scale of clouds from fc m i n = 
20kpc~ 1 (D") to fc min = ^kpc" 1 (D 4 ' ), the feedback 
provided by the jet accelerates the clouds to a velocity a 
factor of two greater, from 600 km s^ 1 to 1200 km s -1 . 
Doubling the cloud sizes from fc m ; n = 20kpc~ 1 to 
fcmin = lOkpc -1 decreases the maximum cloud veloc- 
ities reached in the simulation by a factor of 3, from 
200 km s -1 to 600 km s" 1 . ?7 cr it is therefore more sensi- 
tive to the maximum sizes of clouds than the volume 
filling factor of clouds. Moreover, the scaling between 
fc m in and (v r ^ w ) is nearly linear between k m \ n = 20kpc _1 
and fcmin = 40kpc _1 , and somewhat steeper than linear 
between fc m i n = 10kpc _1 and fc m i n = 20kpc~ 1 . 

The reason for the strong cloud-size dependence and 
linear scaling is that changing the cloud sizes at con- 
stant filling factor changes the rate of ablation relative 
to the total cloud mass without changing the jet plasma 
confinement time. This is because only the amount of 
surface area exposed to ablation relative to the volume 
of a cloud changes. Since fc m i n oc R ~ , where R c is the 
cloud radius, the ratio of surface area to volume of a 
cloud scales linearly with fc m i n . For higher fc m i n , the rate 
of ablation relative to the total mass of clouds increases, 
while the confinement time of the jet plasma does not 
change compared to runs with different fc m i n but iden- 
tical fv- This allows for a far higher fraction of warm 
phase mass to be accelerated to higher velocities, increas- 
ing the maximum value of {v r ^ w ) reached in the run. 

An equivalent statement to the above explanation uses 
the concept of a jet-cloud "interaction depth", Tj C , for a 
given distribution of clouds with varying fc nnn - In anal- 
ogy to optical depth, the effective interaction depth may 
be written tj c = u c -kR\ max -Rbuigc, where n c is the num- 
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Figure 7. Maximum mean radial velocity of clouds, (v r w), 
versus minimum sampling wave number, fc m i n (and correspond- 
ing maximum cloud size, ij c , max ), for runs in the D-series (dia- 
mond points, Pj ot = 10 45 ergs -1 ) and G-series (triangular points, 
Pj ct = 10 44 ergs — 1 ). The runs with filling factor fy = 0.052 are 
marked with black markers, labeled, and connected with lines and 
show the variation of (iv jTO ) with fc m ; n . The grey markers cluster- 
ing around a black marker are runs differing only in filling factor. 
See the electronic edition of the Journal for a color version of this 
figure. 



ber of clouds per unit volume (the number density of 
clouds), and i?buige is the radius of the region in the 
bulge which contains clouds. The clouds may be thought 
of as N scattering centers with a cross-section nRl max 

randomly distributed in the volume (47r/3)i?^ ulgc , and 
the interaction depth may be thought of as a measure 
of the average number of jet-cloud interactions any jet 
stream starting from the origin, including trajectories 
along secondary streams, will experience. This formu- 
lation of an interaction dep th is indeed relevant because, 
as we demonstrate in S5.6 the jet streams carrying en- 
trained hot and warm phase material are directly respon- 
sible for the acceleration of clouds through their ram 
pressure. The total number of clouds in the bulge is 
N = /v-Rbuigc/^max = ^c-Rbulge- Therefore, the num- 
ber density of clouds is n c = fy / R£ max , and the inter- 
action depth is T jc = 7r/v(i? bu lgo/i?c,max) = TTjVfcmin- 

Hence, for a fixed fy, r- ic oc fc mm - 

The linear relation between the ratio of surface area 
to volume of a cloud and fc m i n , or equivalently, the lin- 
ear relation between Tj c and fc m i n , leads to a linear re- 
lation between k mm and (v r , w ), which is seen between 
fcmin = 20kpc~ 1 and fc m i n = 40kpc _1 . It is, how- 
ever, not clear whether one may extrapolate this rela- 
tion to larger cloud complexes with sizes characteristic 
of giant molecular clouds (GMC), say of order several 



100 pc, given that the scaling between k mm — lOkpc 
and fc m in = 20kpc~ 1 is steeper than linear. A possi- 
ble reason for the steepening at larger i? c . max is that the 
larger inter-cloud voids cause a decollimation of the jet 
streams leading to less efficient momentum transfer. The 
scaling may also be affected by resolution limitations to 
capturing the fractal surface of clouds, and by statistical 
variations in the decreased number of jet-cloud interac- 
tions for small fc m i n . It is difficult to predict the feed- 
back efficiency with respect to GMC with scales of order 
100 pc from our simulations because these are generally 
not spherical and the effective interaction cross-section 
depends on orientat ion with respect to the j et stre ams. 
The simulations by Sutherland & Bicknell ( 2007 1 and 
|Gaibler et al.| (2011) show that, it the molecular gas is 
distributed as a large coherent complex in a disc-like ge- 
ometry, the coupling between jet and ISM in terms of 
negative feedback through gas expulsion is weak. Ob- 
servations of some gas-rich radio galaxies indicate that 
the molecular gas is not coupled as strongly into out- 
flows with the jet as the neutral o r ionized material ( Ogle 
eTan2010l IGuillard et al.||2012|), although 4C 12.50~Isa 



prominent exception (jDasyra & Combes 2011 2012). 



The explanations given here also apply to the influence 
of cloud sizes on the mechanical advantage and energy 
transfer efficiency from jet to warm phase in these sys- 
tems. Both the mechanical advantage (Fig|4| and the 
energy transfer efficiency (Fig[5| are significantly greater 
in systems with smaller cloud sizes. 

The sense in which the fractional cloud dispersal rate 
depends on cloud sizes is the same as the dependence of 
the conditions for star formation in a cloud on its size, in 
that, the larger a cloud the more likely it is to collapse 
due to an external pressure trigger. Thus, whether jet 
mediated feedback induces or inhibits star-formation is 
a sensitive function of the statistics of the warm phase 
distribution, in particular its size distribution. 

5.4. The expansion rate of the quasi- spherical bubble 

In this section, we determine the departure of the out- 
flow energetics from that of an energy-driven bubble as 
functions of warm phase parameters. Because our simu- 
lations include radiative cooling and a porous two-phase 
ISM, we expect the energetics of the bubble that sweeps 
up the ISM imparting momentum and energy to the 
clouds to lie between the energy-driven and momentum- 
driven limits. We discuss momentum-driven and energy- 
driven outflows in relation to work in the literature sep- 
arately in |5.5| 

Figure [8^ contains six panels showing the evolution of 
the bubble radius with time for different runs in the G- 
series, D-series, and B-series. We defined the bubble ra- 
dius to be the radius of a hemisphere whose volume is 
equivalent to that swept up by the pressure bubble in the 
simulation. In each panel, the solid black line and the 
dotted black line represent the theoretical, self-similar, 
spherically symmetric evolution of the forward shock ra- 
dius and contact discontinuity, respectively, of an energy- 
driven bub ble (wind) in a unifor m medium in Stage 1, as 
defined by ~ 
|Begelman|| 



That stage represents an adiabatically 
expanding bubble with constant injection power, which 



Weaver et a!] ( |1977| ) (see also §6 |Bicknell fc 

diaT 



in our case is Pj c t. The solutions of the first stage are 
applicable here because radiative losses, although they 
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The spherically equivalent bubble radius as a function of time for simulations with different ISM parameters in the G-series 



(10 44 ergs I ), D-series (10 45 ergs" 1 ), and B -series (10 46 ergs 1 ) of runs. The runs have identical ISM hot phase density (n^ = 1.0cm 3 ) 
but differ in jet power, warm phase filling factor, and maximum cloud sizes. The solid bla ck line and the dotted black line are position o f 
the forward shock and the position of the contact discontinuity in an energy-driven b ubble (Weaver e t al.|1977||Bicknell fc Begelman 1996). 
The dashed black line is the location of the thin shell in a momentum-driven bubble | |Dyso5Jl984| e.g.], a: Bubbles driven by jet feedback 
in lower filling factor environments evolve closer to classical energy-driven bubbles, b: Bubbles driven by jet feedback in halos with larger 
cloud complexes (at constant filling factor) evolve closer to classical energy-driven bubbles, c: Same as a except for fc m ; n = lOkpc - 1 . d: 
Same as a, but for the G-series of runs, e: Same as 6, but for the G-series of runs. /: The dependence of the bubble expansion rate on 
ISM parameters is similar for all jet powers. See the electronic edition of the Journal for a color version of this figure. 



improve th e structural integrity of clouds and their sur 



vival time (Cooper et al. 20091, are energetically unim 



portant in our simulations. The forward shock and dis- 
continuity evolve according to R2 — 0.88(Pj e t/^?i) 1 ^ 5 i 3 ^ 5 
and Rc = 0.86i?2, respectively. The location of the thin 
shell in the momentum-driven limit of a bubble expand- 
ing in a uniform medium of mass density ph is delin- 
eated by the dashed black line, and given by the equation 

-Rsheii = \/3/2 (p/ph) 1 J 4 t 1 / 2 , where p is the momentum 
injection rate ( Dyson|[l984 ) . 

Without focusing on a particular panel in Fig. [8j we 
note that the bubble radius in some runs follows that 
of the theoretical prediction for an energy-driven bub- 
ble closely, while in others the bubble radius initially in- 
creases more slowly than the rate predicted by theory. A 
slight deceleration can even be seen in some runs as the 
bubble is increasingly mass-loaded by warm phase mate- 
rial. In a few runs a return toward the theoretical line is 
visible, with gradients steeper than the theoretical limit 
for an energy-driven wind of given injection power. This 
happens because the medium inside the pressure bubble 
while the jet plasma is confined by clouds is at a higher 
pressure than that of a bubble that is expanding in a ho- 
mogeneous atmosphere with ambient density . During 



the breakout phase of the jet from the region filled with 
clouds, the jet plasma bursts out of the outermost porous 
channels and momentarily fills volumes at a faster rate 
than a bubble that was not impeded and confined for 
some duration by a porous, dense distribution of clouds. 

The first panel (a) shows four runs of differing filling 
factor from the D-series, for which fc min = 20kpc _1 . A 
bubble evolving in a system with larger filling factor ex- 
pands more slowly. As we decrease the filling factor, the 
deviation from an energy-driven bubble become smaller. 
We see the same behaviour for the runs in the G-series 
(panel d). The larger volume of channels available for 
the jet plasma to flood through, and the resulting smaller 
confinement time, is the dominant factor that defines the 
bubble expansion rate. The same trend is visible in the 
second panel for simulations of differing filling factor, for 
which fc min = lOkpc -1 (panel c), although the effect is 
much weaker. This is the result of the confinement times 
and mass loading from hydrodynamic ablation ceasing to 
vary much as the maximum cloud sizes increases. 

The expansion rate profiles for runs with differing max- 
imum cloud sizes, but equal filling factor at fy = 0.052 
is shown in panels (6) and (e) for the D-series and G- 
series, respectively. The expansion rate of the bubble 
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deviates increasingly from the theoretical rate with de- 
creasing maximum cloud sizes. The reason for this is the 
increased mass ablation r ate r elative to the total mass 
of clouds, as described in §5.3| Smaller maximum cloud 
sizes for the same porosity lead to higher mass-loading 
rates and decreased expansion rates of the bubble. 

While feedback efficiencies are mainly sensitive to the 
maximum cloud sizes, the deviation of the bubble expan- 
sion rate from a theoretical energy-driven rate is sensi- 
tive to both the maximum cloud sizes and filling factor. 
Within the parameter range studied here the deviation 
depends more strongly on filling factor than maximum 
cloud size. For systems containing large clouds with 
small filling factors, the bubble evolution approaches that 
of an ideal energy-driven bubble. Thus, it would seem 
that in this limit, these results encourage a sub-grid AGN 
feedback prescription in cosmological models, in which 
energy is injected isotropically into a small region, even 
if the multiphase ISM conditions in the cores of gravita- 
tional potentials are not adequately resolved. However, 
this limit is not the same as that which leads to the most 
efficient cases of negative AGN feedback. The latter is 
attained in the limit of small filling factors and small 
cloud sizes. A distribution of larger clouds, instead, may 
lead to positive feedback, e.g., pressure-triggered star- 
formation. 

Cosmological SPH models commonly invoke negative 
AGN feedback in a single-phase ISM, which essentially 
corresponds to the hot phase in our simulations. The 
heating rate of the hot phase will therefore likely al- 
ways be accurately captured in these models if the fill- 
ing factor is smaller than ~ 0.1. The dispersal of 
warm and cold gas, on the other hand, requires fc m j n > 
20kpc _1 (i? c ,max < 25 pc). Since the bubble expan- 
sion rate does not depend very strongly on the maxi- 
mum size of clouds, however, we assert that negative 
feedback as implemented in cosmological SPH models in 
a single-phase ISM is consistent with negative feedback 
in our simulations with a two-phase ISM if the warm 
phase filling factor is less than ~ 0.1 and the largest 
cloud complexes are smaller than ~ 25 pc. In this regime, 
the embedded warm-phase material is accelerated nearly 
isotropically to the bubble expansion speed within the 
dynamical time of the bubble, ensuring that the nega- 
tive feedback affects both phases, while the bubble re- 
mains approximately energy-driven. This conclusion is 
independent of jet power. 



5.5. Energy- or momentum-driven? 

The theories o f a momentu m-dr iven wind developed 
by |Fabian| (l999[ ) , |King| ( |2003[ ) , and |Murray et al.| ( |2005 l 
naturally predict Mbh oc <t 4 



The theory o f an en ergy - 



driven wind put forward by Silk & Rees ( 1998 1 pre- 
dicts the relation Mbh oc a 5 . These two relations and 
their normalizations are limiting cases for outflow ve- 
locities that can be reached in an outflow powered by 
Ledd(MBH)- In the former case, the outflow loses its in- 
ternal energy through radiative processes (e.g. Inverse 
Compton cooling) and its dynamics is governed solely by 
momentum conservation, while in the latter case energy 
is fully conserved. The difference of 1 in the exponent 
of the relations is not surprising from dimensional ar- 
guments, since energy conservation entails a dependence 



on velocity squared as opposed to a linear dependence 
on velocity associated with momentum-driven flows. 

In FigjSJ the solid lines of constant r\ represent the lim- 
iting slope of a momentum-driven outflow and the dashed 
lines represent the limiting slope for an energy-driven 
outflow. The loci of the maximum values of (v rw ) in our 
simulations with identical filling factor and ISM densi- 
ties and A: m in = 20kpc _1 cluster between log^) = —2 
and log (77) = —4 along narrow strips roughly parallel to 
r\ isolines. The average gradient of the lines connecting 
the loci of (tVytu) appears to lie between 1/4 and 1/5, 
indicating that the outflow is somewhere between mo- 
mentum and energy-driven. 

In a two-phase medium, the determination of whether 
a bubble evolves in the momentum-driven regime or 
energy-driven regime depends on how one compares the 
evolution with that for the case of a smooth, single-phase 
ambient medium. This, in turn, depen ds o n what feed- 
back criteria one is interested in. In 95.41 we saw that 



a bubble evolves in the energy-driven regime as long as 
the warm phase volume filling factors were not larger 
than 0.1. In this regime, the radial heating rate of the 
hot phase (but not necessarily the warm phase) is well 
described by that of an energy-driven bubble. The sup- 
pression of star- format ion in existing clouds is effective 
only if the additional constraint of i? c ,max ^ 25 pc is 
satisfied, because the clouds are then efficiently ablated, 
heated, and dispersed. Only in this regime can the en- 
tire two-phase medium be considered as expanding ap- 
proximately in the energy-driven limit, because the warm 
phase is accelerated to the bubble expansion speed within 
the dynamical time of the bubble. 

One expects energy-driven and momentum-driven out- 
flows to differ in their kinetic power required to achieve 
feedback. For exa mple, consider the predictions of the 



Silk & Rees ( 1998 ) model, which employs the same con- 
dition tor feedback as we have used in our work to derive 
the M-a relation, namely that the outflow velocity ex- 
ceed the host galaxy's bulge velocity dispersion. The 
normalization to the derived M-a relation contains a 
wind efficiency parameter f w — M out Wo Ut /L c dd, where 
u ou t is the outflow velocity, and M out is the mass injec- 
tion rate of the wind. The values of f w can be compared 
to those of rj in this work, because i?kin.u>/-Pjct ~ 0.1 - 
0.4 in our simulations (see Fig. pi). From observational 
estimates of the energetics of AGN outflows, one would 
expect this factor to be of order . 001 - 0.01 (e.g. [McK 
ernan et al.|2007| |Moe et al.||2009[ |Tombesi etaLl|2010 
which is also found in disc-wind simulations (|Kurosawa 
eTal] [2009] |Takeuchi etldT||2lH0l |Ohsuga fc Mincshigc 
201 1| ). This is the level at which cosmological SPH simu- 



LU 



latio ns typically inject energy to model AGN feedback 



Di Matteo et al. 2005 Okamoto 2008 Booth & 



Schaye 2009]) . However, comparing the normalization 



of the observed M-a relation to that derived by Silk 
& Rees, one obtains f w ~ 7 x 10 _6 / gas , where / gas is 
the gas fraction in the dark matter halo, indicating that 
in the spherically symmetric energy-driven regime, low 
Eddington ratio outflows arc sufficient to significantly 
disperse gas. This is also the result found in some semi- 



analytic models (e.g. iCroton et al. 20061) where heating 



by AGN with small Eddington ratios suffices to suppress 
star-formation in massive galaxies and offset cooling in 
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clusters. A possible reconciliation for these inconsisten- 
cies is that the observed outflows are momentum-driven. 
For example, a momentum-driven wind maintained by 
the photon momentum flux of an Eddington limited ac- 
cretion flow (Ledd/ C ) requires f w — v out /c ~ few x 1CP 3 
to accelerate a spherically symmetric shell of swept up 
ISM to a velocity v = a. 

As mentioned however, one must be careful about the 
types of feedback cri teria one is comparing in the vari- 
ous studies (see also { 6.2 1, and how simulations partition 
the injected energy into thermal or kinetic energy in sub- 
grid fe e dback p r escrip ti ons. The theories by |Silk & Rees| 



faylSn 

fll998| ), |Fabian| ( |1999[ ), |King| ( |2003[ ), an d [Murray et al 



( 2U(J5p aim to explain the M-a relation. To inhibit star- 
formation within a galaxy it suffices to merely heat or 
ionize the dense gas. To offset cold gas accretion and 
avoid late-time star-formatio n, more powerful ou tflows 
are required to heat the IGM. Silk fc Nusser| ( [2010[ ) argue 
that, Eddington-limited AGIS do not provide enough en- 
ergy during their lifetime to generate momentum-driven 
or energy-driven winds that can unbind the gas from the 
galaxy potential. AGN jets arc sometimes favoured in 
these cases. In cosmological SPH simulations, AGN feed- 
back is usually implemented as thermal energy injection 
into particles near the core of galaxies, which effectively 
results in an energy-driven bubble with maximal mechan- 
ical advantage that heats the ISM and efficiently inhibits 
local star formation and cluster-scale cooling flows. The 
feedback requires relatively high injection rates of ther- 
mal energy, compared to a wind that is dominated by 
kinetic energy and can drive the bubble through ram 
pressure as well thermal pressure (lOstriker et al. 2010k. 

In this work, we found that the minimum value of r\ 
required by a jet to disperse the warm phase to the veloc- 
ity dispersion implied by the M-a relation is r\ > 10~ 4 
and depends on the ISM density, filling factor, and cloud 
sizes. The regime of fy > 0.1 and i? c ,max < 25 pc comes 
closest to an energy-driven bubble of the entire two-phase 
medium, and, ?7 cr ;t is accordingly small (< 10~ 4 ). This 
is not surprising since the limit fc m j n — > oo is essentially a 
single-phase medium, and the surface area for ram pres- 
sure and thermal pressure to work on (and consequently 
the mechanical advantage) is maximised. Even though 
the bubble expansion rate depends quite sensitively on 
filling factor, ?7 cr it does not. For i? c ,max ^ 25 pc the 
clouds are not strongly ablated and accelerated by the 
bubble. The bubble containing mainly hot phase gas ex- 
pands in a nearly energy-driven manner, but 77 cr ;t is large. 
These results for 77 cr j t , therefore, demonstrate that both 
filling factor and the size distribution of clouds need to 
be considered when assessing the efficiency (and type) 
of feedback, because these factors influence the degree 
to which the warm phase is incorporated in the outflow. 
Conversely, whether feedback is efficient or not cannot 
be uniquely determined by assessing whether an outflow 
is energy- or momentum-driven. 

In our simulations the bubble evolves near the energy- 
driven limit unless fv ^ 0.1. Other indications support- 
ing this approximation are: 1) the mechanical advantage 
(Fig. [4]) and the efficiency of energy transfer from jet to 
the warm phase (Fig. [5]) are high at all times; 2) For con- 
stant /c m ; n , the feedback efficiency scales roughly with 
Pjet/nh (Fig. [6j), a characteristic parameter of energy- 
driven outflows. In the following section we investigate in 



detail how the warm-phase material is accelerated nearly 
isotropically to the bubble expansion speed within the 
dynamical time of the bubble, under the assumption that 
the bubble evolves in the energy-driven regime. 

5.6. Cloud acceleration through ram-pressure driving 

We have investigated in more detail the physical mech- 
anisms that accelerate the clouds. As contained in the 
fluid equations (Eqns.[T]), gradients in both ram pressure 
and thermal pressure contribute to momentum transfer 
from the jet plasma to the warm phase. The mechanical 
advantages measured in WB11 and here very high (s> 1) 
for all simulations (see Fig. [4]). 

We first show that the simple picture that the expan- 
sion velocity of the jet-driven bubble exerting a ram pres- 
sure on the clouds is responsible for the cloud velocity 
does not work. Careful inspection of the flow in the sim- 
ulations, instead, shows a combination of other effects 
combining to provide a high mechanical advantage. 

The analytic expression for a jet-blown bubble in the 
energy-conserving limit is not a bad approximation in 
this context as has been shown by our simulations. The 
radius R}, of a bubble blown by a jet with power Pj et into 
a medium with ambient density p a is given as a function 
of time t by 

Rb = At 3 / 5 , (7) 

where 



.4 



/125R 



V3847rp a 



1/5 



(8) 



The driving of a spherical cloud of mass m c = Air/3R 3 



and radius R c , to a velocity v c 



via the ram pres- 



sure of the expanding bubble, which has rest mass den- 
sity pb and expansion velocity i> eX p = (3/5)At~ 2 / 5 is de- 
scribed by the equation of motion: 



' dt 



(9) 



where Crj is the drag coefficient and ~kR\ is the cross- 
sectional area of the cloud. The acceleration of a cloud 
thus, 



is 



dv c 
~dl 



3C D p b v, 



cxp 



(10) 



4 p c Rc 

If we assume in this initial model that the density of 



the bubble is determined the jet mass flux M\ ct , then 



d 
dt 



4tt 

Y 



PbRb 



and, integrating, the bubble density is 

Pb = ;'' ,1 :i M,:i ■ ' 



(11) 



(12) 



We can use equation ( 10 1 to calculate an acceleration 
timescale for the cloud from 



"^GXp 3 

= — Ox? . , _ 

tucC 4 V Pc / R' 



2 

cxp 



(13) 



and this implies that the acceleration timescale for clouds 




(14) 
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The acceleration timescale for clouds compared to the 
evolution time for the bubble is: 



^acc 



807T i 

~W Cd 



RrAH 1 '* . 



(15) 



jet 



In order to evaluate £ acc we need to determine the jet 
mass flux. This can be determined from the jet power as 
follows. For a jet with relativistic enthalpy, w, velocity 
c/3 and cross-sectional area Aj et the jet power is given by 



P 



jot 



T 2 wc/3A jct - Mc 



(16) 



where the mass flux, M, for a jet with proper rest-mass 
density is given by 

M = TpcpA^ t . (17) 

The energy flux equation has two ter ms. The first term 
is the conventional relativistic form (Landau & Lifshitz 
19871 which incorporates the energy flux due to the rest 



mass energy plus the normal kinetic and internal energy 
terms. These contributions all originate from the rela- 
tivistic enthalpy w. The second term subtracts off the 
rest mass energy flux since, in this context this is not 
useful energy. (Nuclear reactions are not involved.) 

The ratio of the dynamic variables appearing in the 
energy and rest-energy fluxes is 



T 2 vc(3 
Tpc 3 /3 



w 



pes 



Hence the energy flux is given by 



Pjet — 



1 I M jet c 



(18) 



(19) 



We put 



w = pc 2 + (e + p) , 



(20) 

where e is the internal energy density and p is the pres- 
sure. This gives 



Pi 



jet 



(r 



DMj ct C 2 



where \ = P° 2 /( e 
expressed ; " ' ' 



1 



1 



(21) 



-p). Therefore, the mass flux may be 
in terms of the energy flux by: 



Pjet 



r - 1 c 2 



1 



r-ix 



(22) 



For typical values, e.g., T = 10, Pj 0t = 10 45 ergs~ 
X = 1.6, Cd = 1, p c = 1000, p a = 1, and R c — 25 pc, it is 
obvious that this model does not work because t acc /t ^> 
1. 

Closer inspection of the ram pressure vectors (p|u|u), 
and thermal pressure gradients in the flow field in Fig 
[9] reveals a combination of effects that accelerate clouds 
faster than on a bubble evolution timescale. In the di- 
verted secondary jet flow channels we do not see the large 
thermal pressure gradients maintained along the primary 
jet axis when the jet head encounters a cloud. It is evi- 
dent from inspection of the density maps in Figs. [T] and [9] 
that the mean density in the plasma flow that is flooding 
through the channels between the clouds is much larger 
than the mean density of a purely jet-blown bubble. Here 



the transfer of momentum is primarily maintained by a 
ram pressure that is comparable or somewhat greater to 
that of the primary jet flow by virtue of turbulent mix- 
ing (entrainment) of the shocked hot-phase gas with the 
jet plasma, and hydrodynamic ablation of cloud mate- 
rial into the engulfing flow. In channels with high mass- 
loading, the ram pressure even exceeds that of the pri- 
mary jet. 

Another quan tita tive discrepancy with the theory con- 
tained in Eqn. (15) is that the flow velocities providing 



the ram pressure are not the expansion velocity of a fully 



thermalized bubble, 



u GXp 1 



but the velocity of partially 



thermalized channel flow, v c h ~ 10 kms . 

Thus, two modifications need to be made to the above 
theory: 1) The density of the bubble in which the clouds 
are embedded is not solely determined by the mass in- 
jection of the jet, but the mass injection needs to be 
enhanced by the entrainment rate of hot phase material 
and the ablation rate of warm phase material; 2) The 
channel flow velocity that carries the momentum and 
provides the ram pressure at channel-cloud interfaces is 
not the expansion speed of the bubble but a much higher 
speed of only partially thermalized material. With these 
additi ons to the theory, we derive modified versions of 
Eqn. ([15). 

In regard to point 2) abov e, we write v c h instead of 



on the RHS of Eqn. ( 13 ) 



-C 



D 



R C 



(23) 



If the channel flow in the bubble is entraining hot- 
phase material and mass-loaded by hydrostatic ablation 
of clouds, the mass injection into the bubble is aug- 
mented to 

M tot = Mjet + Mentr + M abl , (24) 

where the entrainment rate is approximately equal to the 
rate at which matter is swept up by the spherical bubble, 



Mentr = 47r/ ontr p a P 2 W oxp , 



(25) 

given the fraction of entrained material, / on tr- For the 
case of purely hydrostatic ablatiorj^J the ablation rate of 
one cloud with internal isothermal sound speed c r embed- 
ded in a cha nnel flow with Mach number M. c h (Hartquist 
et al.||1986[ ) is: 



amin (l,M*Jf \ (m c c c ) 2/:i {p b v ch ) 



1/3 



(26) 



and the total mass injection rate into the bubble ablation 
rate is given by 



M a bi = rh, 



x —R 3 f„k 3 



(27) 



The constant a is of order unity for spherical clouds. 

The treatment that follows is not entirely self- 
consistent because we do not determine the bubble ex- 
pansion rate, Rb(t), self-consistently under the modified 
conditions, but continue instead with the assumption 
that the outflow remains close to an energy-driven bubble 

4 The ablation in this model is driven by pressure differences 
developing around the cloud surface, with ablation rates highest 
perpendicular to the flow direction. 
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Figure 9. Density and pressure maps of the mid-plane region near the jet base during three stages of jet-cloud interactions in run D'. The 
upper row shows a color map of the density. In the middle and lower rows, ram pressure vectors (pv|-u|) and thermal pressure contours (p) 
are superimposed with a color scaling in units of \i cm - 3 c 2 (fj, is the mean mass per particle and c is the speed of light) on a greyscale density 
map. The left, center, and right columns show the data at various times at which different effects dominate the acceleration mechanism of 
clouds: Left: the jet head is strongly interacting with the clouds in its path; Middle: hot-phase entraining jet streams carry ram pressure 
to clouds embedded in the bubble; Right: jet streams carrying entrained hot-phase and ablated cloud material dominate the channel flow. 
Compared to the ram pressure, the thermal pressure is relatively uniform inside the bubble. See the electronic edition of the Journal for a 
color version of this figure. 

that follows Eqn. fcf\. In 1 5.4 we saw that this approxi- Let us look at the case where M e ntr ^ Mj et and M Gntr 3> 
mation is reasonable. The bubble density is obtained by 
integrating: 



d 

dt 



47T 



Pb Ri = M je 



M, 



a hi 



(28) 



Ma.hU which holds early in the evolution of the bubble. 
We retain only the second term on the RHS in Eqn. ( 29 1 
and substitute that expression for pb into Eqn. ( 23 ) to 
obtain: 



t - V 1 Pc v * 

3 JcntvPa V 



(30) 



ch 



where M jet , M, 



on t r , and Mabi are defined by the equations 
above. This gives 



The acceleration timescale as a fraction of the dynamical 
time is then 



Pb = ^A- 3 M ict r^ 5 + f cntlPa + m c f v k 3 min ^t . (29) 



t. 



acc 

T 



Rc 



fentrPa 



At~ 7 / 5 . 



(31) 
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Equation (31) shows that clouds embedded in a heav- 
ily entrained spherical bubble will experience more effi- 
cient acce leration with time. Inserting typical values into 
Eqn. pi] ), in particular, / entr = 1, and v ch = w c h,entr = 
10 5 kms _1 , we obtain t acc /t ~ 0.078 at t = 100 kyr. 

At t = 100 kyr, the bubble density is, however, domi- 
nated by the mass loading from cloud ablation. We there- 
fore turn to the limit M a bi ^ Mj et and Af a bi ^> M cntr . 
Retai ning only the third term and using Eqn. ( [26] ) , 
Eqn. (29) becomes 



p b = a min (^1,X^ 3 ) (m c c c ) 2/3 (pbV c h) 1/3 fvk^—t 



47T 

T 



= ^RlPcCc min (l, M 2 ch ) v\£ 



" 1 Yi a ^ yk 



3/2 

3 \ ,3/2 
min ' 



(32) 



We may estimate the channel speed from the (local) mass 
continuity condition that the mass ablation rate reach 
the channel mass flux through an area irR 2 , m c ,m a x = 
Pf>Vch,abi7r#c, in which case, using p b = TO Ci /y&4 in (5/14)t 
(Eqn. 29} , 

14 -i 
v*, a u= 5nR2Jvk3 t . (33) 



Combining Eqns. (|23|) , 132I) , and (33), we obtain 



Pb - 



(ttR 2 ) 1/2 R c p c c c (\ 



3/2 



XM 2 ch )f 



rh 6 f 



10 
21 

x min 



xmax(l,7W ch 2 ) (f v k 3 mi 



TTRt 



,3/2 



-3/2 



At 



-2/5 



(34) 



(35) 



Note that the expression for i acc /t (Eqn. 35) only de- 
pends on the isothermal sound speed in the cloud, rather 
than directly on the density. In our simulations we ob- 
serve that the channel flow dominated by entrained ma- 
terial moves with internal mach number, A4 ca ~ 1, while 
the channel flow dominated by ablated cloud material 
moves with internal Mach number Ai ca ~ 3. Taking 
c c = y/kT c /n with T c = 100 K, f v = 1, fi m = 0.6165, 
and /c m in = 20kpc~ 1 , we find t &cc /t = 15.7 for a = 1 and 
iacc/* = 0.49 for a = 10 at * = 100 kyr. 

In Fig. [10] we show the curves for the expressions for 
Pb and iaccA obtained above. Unless otherwise men- 
tioned we used the set of typical parameters mentioned 
above. In the left panel we show the predicted mean 
bubble density for the case for which Mj e t ^> M cntr and 
Mj et > M ab i (Eqn. 12 with Eqn. 22 ) with a dotted blue 
line, the case for which M ontr ^> Mj et and M entr 3> M a bi, 
that is, pi, — f entr Paj with a dash-dotted green line, and 
the case for which M a bl 2> Mentr and M a bi 3> Mjet 



(Eqn. 34 ) with dashed orange lines. We see immediately 
that the maximum contribution to the total density by 
the jet plasma alone is much smaller compared to that 
of the other mass injection mechanisms. Since hot-phase 
entrainment is effective from t = 0, jet plasma mass in- 
jection never dominates throughout the evolution of the 
bubble. Instead, the entrained material mixes well with 



the jet plasma, and the average bubble density is close 
to the ambient density. 

In time, cloud ablation becomes important in con- 
tributing to the total density. The bubble age at which 
this happens depends on p c , fc m ; n , R c , fy, and a, and 
the t herm odynamic warm phase parameters. We plot 
Eqn. (34) for a = 1, a = 10, a = 20, and for the case 
for which fc min = 40kpc~ , R c = 12.5 pc, a = 10. Su- 
perimposed are points from Simulation D', subdivided 
into the phases before, during, and after jet breakout. 
The error bars denote the upper and lower limits in the 
estimate of the mean bubble density in our simulations 
obtained by choosing different cutoffs for the tracer vari- 
ables. The curves tracing the summed contr ibut ions to 
the total density are obtained by solving Eqn. ( 29 ), which 



is an implicit equation in pb and t, if all three terms on 
the RHS arc included, with a standard root finding algo- 
rithm. 

We find that cloud ablation is quite efficient in our 
simulations, requiring a > 10 to match the mean bubble 
densities seen in our simulations. Hartquist et al. ( 1986 1 
introduced a as a constant near unity in the ideal case m 
which embedded clouds are spherical. Spatially fractal 
clouds have an (undefinably) larger surface area exposed 
to ablation, which leads to a larger inferred value for a 
from our simulations. The fact that the analytic curve 
for the case for which fc m i n = 40kpc~\ R c = 12.5 pc, 
a = 10 is similar to the curve for which fc m j n = 20 kpc -1 , 
R c = 25 pc, a — 10 lends support to the theory that the 
total surface area exposed to ablation by the channel flow 
is the main parameter governing the global bubble evo- 
lution, since a cloud distribution with fc m i n = 40kpc~ 1 , 
R c = 12.5 pc contains the same mass as a cloud distri- 
bution with fc m ; n = 20kpc~ 1 , R c = 25 pc but four times 
the surfa ce a rea for ablation. 

In Fig. 10 b) we plot the acceleration timescale, t acc /t 
against t, for the same limiting cases as in Fig. 10 a). To 
evaluate the acceleration timescale when all mass injec- 
tion terms are included, we use the solution to Eqn. (29 1 
and Eqn. (33) in Eqn. (23). Because of the strong ae- 
pendence ot t he acceleration timescale on channel veloc- 



ity (Eqn. 23), we take into account that Eqn. (33) is 
only valid once cloud ablation becomes important T5y in- 
troducing an exponential turnover from v ca = i> c h,entr to 
^ch = v c h abi at t — o ver a timescale defined by equating 
the RHS of Eqn. ([35]) with f entI p a - 



The acceleration timescale predicted by the model of 
clouds embedded in an expanding bubble (dotted lines) 
is too large to explain the cloud velocities seen in our sim- 
ulations. The acceleration timescales predicted by ram- 
pressure driven, entrained and mass-loaded channel flows 
(solid lines), on the other hand, fall below t acc /t = 1 after 
~ 10 kyr for all a. In the early phases of the bubble evolu- 
tion, the acceleration of clouds quickly becomes efficient 
as the hot-phase entraining jet plasma provides a high 
ram pressure to the clouds. The acceleration of clouds 
becomes somewhat less efficient when mass-loading from 
cloud ablation becomes important, as this reduces the 
channel flow velocity. 

The analytic estimates above only represent the 
global mean bubble density and mean cloud acceleration 
timescales in a spherically symmetric bubble, whereas, 
in reality, there exist radial and angular variations of the 
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Figure 10. a: The mean density of an energy-driven bubble as a function of time. Individual contributions to the density are shown in 
separate broken colored lines, as indicated in the legend. The combined contribution is shown in solid black lines for different values of a 
and fc m i n . Superimposed are data points from Simulation D'. The error bars denote estimate limits using different cutoff values for the 
tracer variable, b: The ratio of cloud acceleration timescale to bubble age, versus time. The acceleration timescales as predicted for the 
cases that individual contributions dominate the bubble density and channel speed are indicated in broken colored lines. Note that the 
curves for a = 10 and for fc m i n = 40kpc — 1 , a = 10 for the cloud-ablation dominated case (dashed lines) overlap. The solid lines trace 
the (approximate) combined mean acceleration timescale, taking into account the transition from an entrainment-dominated bubble to a 
cloud-ablation dominated bubble. See the electronic edition of the Journal for a color version of this figure. 



mean bubble density and channel speed, while clouds 
are driven outward. The relative contributions of jet 
plasma, hot-phase entrained channel flow, and warm- 
phase mass-loaded channel flow are, however, reasonably 
well explained by the analytic description when one al- 
lows for a high value of a due to the greater surface area 
of fractal clouds exposed to the ablating channel flow. 
Improvements to the analytic description may include 
a self-consistent evaluation of i?&(t) and a finite spatial 
cloud distribution and finite cloud masses available for 
mass loading the channel flow, but such a treatment is 
beyond the scope of this paper. 

We conclude that it is the ram pressure of mass-loaded 
channel flow, resulting from entrainment of the ambient 
medium and ablation of clouds, which provides the pres- 
sure gradients at cloud interfaces that transfer momen- 
tum and result in sustained high mechanical advantage 
and efficient bulk cloud acceleration within the dynami- 
cal time of the bubble. This acceleration mechanism of 
clouds found in our simulations may be interpreted as a 
variation of the t wo-st age feedback model proposed by 
Hopkins fc Elyia| p010| . 

Hopkins &: Elvis proposed a radiative mode of quasar 
feedback, in which a radiation-driven quasar wind of 
the hot phase engulfs the warm and cold phase clouds 
in a galaxy and transfers energy and momentum to 
the clouds, accelerating and disrupting them through 



sustained pressure gradients and dynamical instabili- 
ties. The enlarged surface area of ablated and expanded 
clouds provide a larger cross-section for photoionization 
and radiation pressure from the quasar. Overall, this 
means that the efficiency of negative feedback is higher 
than that in the case of direct irradiation of clouds by 
the central AGN. The scenario the authors describe is 
very similar to the flow evolution and acceleration mech- 
anisms we identified in our simulations. The jet plasma 
flow with its entrained hot phase material is analogous to 
the hot-phase quasar wind, and, as the flow engulfs and 
ablates clouds, sustained pressure gradients drive out the 
cloud material. In the case of our simulations, the pres- 
sure gradients are maintained by ram pressure, which 
may also be the case for quasar winds, in addition to 
radiation pressure. 

6. DISCUSSION 

6.1. The feedback efficiencies in radio galaxies with 
detected outflows 

There are a growing number of observational studies of 
radio galaxies in which outflows of cold neutral or warm 
ionized material are detected. We have compiled a set 
of 27 radio galaxies from the samples studied by |Holt 



(2008), Nesvadba et al.| ( |2006l |2008 2010), 



uillard ct al. (20 1 21, a nd the individual 



Torresi 



case 



Stockton et al. 1 5 



We also use the 
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Table 3 

Radio galaxies with outflows 



Radio galaxy 


log Pjet (a) 


v out (ionised) 


(b) iw (HI) M 




M BH (C) 


Ref( f ) 


v (s) 




(erg) 


( km s ) 


(kms" 1 ) 




(Mq) 






|Holt et al. ( 2008 ) 


3C 213.1 


45.lt 


142 




0.19 


9.1 


W09 


0.0074 


3C 268.3 


45.9+ 


760 




0.37 


7.8 


W09 


0.93 


3C 277.1 


45.7+ 


79 




0.320 


7.6 


W09 


0.89 


4C 32.44 


46.5* 


852 


128 


0.368 








PKS 1345+12. . 


44.6+ 


1980 


400 


0.122 


7.8 


W09 


0.052 


3C 303.1 


45.5+ 


438 




0.27 


8.1 


W09 


0.11 


PKS 1549-79 . . . 


45.3+ 


679 


79 


0.152 


8.0 


W09 


0.15 


PKS 1814-63 . . . 


45.6+ 


162 


24* 


0.065 


8.7 


Mil 


0.058 


PKS 1934-63 . . . 


45.6+ 


93 




0.182 








PKS 2135-20. . . 


46.8* 


157 




0.636 








PKS 2314+03. . 


45.8+ 


497 


350 


0.220 


8.5 


W09 


0.15 



MRC 1138-262 1 
MRC 0316-257 2 
MRC 0406- 244 2 
TXS 0828+193 2 
3C 326 3 



Nesvadba et al. 



(2006 



2008 



2010) 



1,2,3 



46.2 
46.0 
46.2 
46.4 
45 



800 
670 
960 
800 
1800** 



2.16 
3.13 
2.44 
2.57 
0.09 



8.7 
8.3 
8.5 
8.7 
8.6 



N06 
S07 
S07§ 
S07§ 
HR04 



0.16 
0.40 
0.40 
0.40 
0.020 



10000 ulv 
600 

43769 DW 

1000 



Torrcsi et al. (2012) 



3C 445 .. . 
3C 390.3 . 
3C 390.3 . 
3C 382 .. . 



44.9 
45.1 
45.1 
44.8 



0.0562 
0.0561 
0.0561 
0.0579 



8.3 
8.6 
8.6 
9.1 



M04 
WU02 
WU02 

M04 



0.028 
0.029 
0.029 
0.0044 



Guillard et al 



3U 236 

3C 293 

3C 305 

3C 459 

4C 12.50 

IC 5063 

OQ 208 

PKS 1549-79 . 



(2012) 



45.9 
45.7 
45.8 
46.4 
45.4 
45.3 
43.7 
45.9 



507 
494 

372 
812 



906 



75LT 

500 

250 

300 

600 

350 

600 

250 



0.10 

0.045 

0.042 

0.23 

0.12 

0.011 

0.077 

0.15 



8.0 
7.9 
8.5 
7.8 
7,1 

8.0 



W09 
WU02 
W09 
W09 
V10 

W09 



0.40 
0.63 
0.63 
0.31 
0.63 

0.63 



3C 48. 



46.5+ 



491 



Stockton et al. (2007) 



0.369 



W09 



0.39 



Jet power. + Values from |Wulf2009| Table 1) *Jet powers neither given in reference or in Wu| l |2009| were computed using the 1.4 GHz 
flux and the scaling relation by Cavagnol o et al.| ( |2010 i. 

( b ' Outflow velocity of ionized gas. "Terminal velocity estimated from Na D absorption blue shift of 350 kms -1 . DW/ Disk wind. 



Revised by Morganti et al. (2011 1 



( c ) Outflow velocity of neutral gas, seen in H I absorption 
< d ) Redshift. 
SMBH mass. 

( f ) Reference or calculation method for SMB H mass. W09: SMBH mass taken from 
from stellar bulge mass, M*, given in |Seymour et aL] |2007 | S07) or|Haring 

Mrh = 0.0015M*. §M t is only an upper limit; M04: Marchesini et al.J |2004|; WU02: [Woo fc Urryj l 2002| ; V10: V.- U+ van H a] 



Wu (2009 



Table 1); S07, HR04: SMBH mass estimated 
HKp4), respective ly, usi ng the Magorrian relation 

~'j2010); 



Mil: |Morganti etHlpbTT) ; N06: 

using the stellar mass estimated by |lNeBvaclDa et al.| l |2O0b"| l and tne Magorrian relaT 
( g ) Jet Eddington ratio, r\ = Pjet/i'edd estimated using Mbh- 



Mr. 



U.U015A1* 



sample of Lehnert et al. (2011 data obtained through 
private communication) in the tollowing comparison with 
our simulation. The galaxies are listed in Table [3] whose 
columns include the jet power, the outflow velocities in 
either neutral or ionized gas (or both), the SMBH mass, 
and the value of r\ = Pj e t/£ e dd estimated from the SMBH 
mass. 

When jet powers were not given in the references 
above, they were either found in the table of pr operties 
(Table 1) of the sample compiled by Wu (2009) or cal- 
culated from the 1.4 GHz flux available from the NASA 



locities are usually accurate to a few tens of percent, the 
velocities only represent a particular phase of gas mov- 
ing at the given bulk velocity. The accuracy of the black 
hole mass estimates depend on the method used but, for 
this sample, they are reliable to within a factor of a few. 
Given these uncertainties, it is not possible to arrive at 
strong conclusions, but the method of comparison itself 
is of some interest. 



Extragalactic Datab ase using the scaling relation by Cav- 
agnolo et al] fl2010D 



We note that the none of the estimates for any of these 
quantities are very accurate. The jet power is typically 
only an order of magnitude estimate. While outflow ve- 



Holt et al. (20081) studied the optical emission line gas 
kinematics m 14 compact radio sources and found strong 
evidence for fast outflows in 11 cases. 8 of the 11 cases 
are also known to have blueshifted Hi absorption. We 
list both neutral and ionized outflow velocities in Ta- 
ble |U The two fastest outflows at ~ 2000 km s -1 and 
~ 850 kms -1 were found in GPS sources, and the au- 
thors report a general trend in their sample that the 
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larger outflow velocities were seen in more compact ob- 
je cts, although orientation also p lays a role. 



ISM. 



Nesvadba et al.| ( 2006 2008| ) observed four HzRGs 
(z ~ 2 - 3) with the SliNFONI integral spectroscopy unit 
on the VLT and detected red- and blueshifted O in emis- 
sion aligned with the jet axis. These galaxies contain 
large ionized gas masses (few 1O 1O M0) comparable to 
the mass of the neutral and molecular gas. The inferred 
outflow speeds of ionized gas were ~ 600 - lOOOkms -1 
and energy coupling e fficiencies between the j et and out- 
flow were of order 0.1. INesvadba et al. (2010) posit that 
star formation in the comparatively nearby radio galaxy 
3C 326 is maintained low by the energy input of the ra- 
dio jet, which drives out a fast outflow of ionized gas 
and keeps the neutral and molecular gas warm, in what 
is sometimes termed the "maintenance phase" of AGN 
feedback. HzRGs, on the other hand, may be in an "es- 
tablishment phase" of galaxy formation, in which stellar 
populations are born in starbursts that are then abruptly 
s hut down by the AG N jet. 



Torresi et al. (2012) describe the properties of the only 
three broad line radio galaxies in which outflows of typi- 
cally 100 to 1000 kms" 1 have been detected through ab- 
sorption lines in the soft X-ray spec trum. Thes e AGN 



are said to contain "wa rm absorbers" ( Blustin et al. 2005 - 
McKernan et al.|2007 ) . They find two sources harboring 



outflows originating in the torus or narrow line region 
rather than being associated with a disc wind. Their 
data extend the tentative positive correlation for Type 1 
QSOs between radio loudness and warm absorber mass 
outflow rate. We include these two sources in our sample 
in consideration of the possibility that the acceleration 
of the outflows is driven by jet-ISM interactions. 

Warm absorbers are observationally distinct from the 
ultra-fast outflows (UFOs) detect ed in radio quiet AGN 



at high incidence (at least 35%, Tombesi et al. 120101 
through highly ionized and blueshift ed Fe K-shell absorp- 
tion lines in the hard X-ray band (Cappi 2006). UFOs 
are thought to be mildly relativistic disc winds (~ 0.01c 
to 0.1c, several 10 3 kms _1 to 10 4 kms _1 ) with mass out- 
flow rates comparable to the accretion rate, and outflow 
kinetic luminos ities a significant fract ion of the bolomet- 
ric luminosity (Tombesi et al. 2012). The discovery of 
UFOs in radio-loud AGN jTornbesi et all 120101 12011 1 



blurs the distinction between disc winds and jets and 
motivates research into unified theories of disc-wind-jet 
systems. 3C 390.3, for example, harbors a UFO with out- 
flow velocity (0.146 ± 0.004)c, in addition to the warm 
absorber. Studies of jet/wind-disc in teractions have fo- 
cused more on galac tic microquasars ( |Fender et al.|2004 



Neilsen &: Lee|2009[), often regarded as scaled down ver 
sions of AGN jjVlirabel fe Rodriguez||1999[), and many 



simulation models are a pplicable to~both microquasars 
and A GN ( |Takeuchi et al.||2010| |Ohsuga fe Mineshige 
20111. UFOs have kinetic luminosities comparable to 
those of AGN jets, so they may have a significant im- 
pact on the galaxy-scale ISM of the host. They tend 
to be much less collimated and somewhat slower than 
jets at the same radius, but the physics of the interac- 
tions with the ISM may be similar to that presented in 
this study. This needs to be verified with commensurate 
simulations of ultra-fast disc winds expanding from sub- 
parsec to galactic scales in an inhomogeneous two-phase 



Guillard et al. (2012) have detected outflows in ionized 
gas through broad blueshifted [Neil] emission lines in 
Spitzer IRS observations of a sample of radio galaxies, 
for which previous studies found neutral outflows in Hi 
absorption. We list both neutral and ionized outflow 
velocities in Table [3j Five of the sources have highly 
blueshifted wings on the [Ne il] line that match with the 
b lueshifted broad H I a bsorption. 

|Lehnert et al. (2011 1 fitted the NaD absorption feature 
of 691 SUSS sources with extended radio morphologies, 
redshift z < 0.2, and 1.4 GHz fluxes greater than 40mJy, 
finding modestly blueshifted but highly broadened ab- 
sorption excesses for about half the sources. The au- 
thors deduce the presence of outflows with a distribution 
of terminal velocities between 150 kms -1 - 1000 kms -1 , 
mean mass and energy outflow rates of lOMQyr -1 and 
1 42 ergs~ 1 , respectively . 



Stockton et al. ( 2007 ) mapped the O in emission near 
the central region of 3C 48, a powerful z « 0.369 
CSS source and ULIRG born from a major merger 



( Scharwachter et al. 2004 ), using the GMOS integral field 
unit on Gemini North. The O III emission is blueshifted 
by ~ 500 kms -1 and offset by ~ lkpc northward of the 
quasar, along the jet axis. It is, therefore, distinct from 
the AGN narr ow line region. The st ellar age estimate 



of the host by Stockton et al. (20071 disfavors the hy 



pothesis that the current jet activity triggered the star- 
burst, but the energetics and the alignment with the jet 
of the outflowing gas support the view that a substan- 
tial amount of material is driven out by the jet as it is 
breaking out of the dense central region. 

Two sour ces feature in bo t h the samples of 
(20081) and iGuillard et al 



ISO] and PKS 1549-79. 



(2012), PKS 



Holt et al. 
1345+12 (4C 
he estimates for the outflow 
velocities are different in each study so we list them here 
separately. In addition to the outflow observed in the 
neutral phase through Hi abso rption (jM organti et al 



20041 and in the ionized phase ( Spoon et al.||20d5 ), the 
TJCTRG and GPS source 4C 12 50 exhibits outflows of 
~ 600 - 1000 km s -1 in warm and cold molecular gas 



( |Dasyra fc Combes 2011, 2012). The observations for 
the molecular outflows are not included in the table but 
the data suggest that 4C 12.50 may be a remarkable case 
of a radio AGN in which the all phases of the ISM are 
strongly coupled to the relativistic jet. For sources for 
which a black hole mass was found in the literature we 
calculated the value of r) = i^et/Asdd- 

We present the data from the compiled observations 
together with those from our simulations on the plane 
defined by outflow velocity and jet power in Fig. [IT] The 



data by Lehnert et al. ( 201 lh are plotted as filled con 



tours in probability density, d-ZV/dPjetdiw- The orange 
contours show the probability densi ty if the scaling be- 
tween 1.4 GHz fl ux and jet power by Birzan et al. (2004) 
is used (see also |Best et al.|[2006 ), and the re d contours 
show the p robability density if the scaling by Cavagnolo 



et al. ([2010[ ) is used. The data of the galaxies compiled in 
Table [3] arc shown with different symbols as indicated in 
the legend. Black markers denote outflow velocity mea- 
surements for ionized gas, while magenta markers indi- 
cate velocities for outflows of neutral gas. The points 
are labeled with the value for r/ = Pj e t/£ e dd; for the 
cases for which the mass of the SMBH was found in 
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Figure 11. Simulation results and data from observations com- 
piled in Table [3] on the plane defined by the outflow speed and 
jet power. Simulation results are shown in blue and green points 
corresponding to ISM hot phase densities «/, = 0.1cm -3 and 
1.0cm -3 , respectively. Triangle, round, and square markers de- 
note fc m i n = lOkpc -1 , 20kpc~ 1 , and 40kpc~ 1 , respectively, while 
the shading of the marker from filled, through dark and light grey, 
to white denote filling factors of 0.42, 0.13, 0.052, and 0.027, re- 
spectively. The samples studied by the various authors are marked 
with different symbols and contours as shown in the legend. The 
black symbols mark the measured outflow speeds of ionized mate- 
rial, and the magenta symbols mark the measured outflow speeds 
of neutral material. Subscripts denote the value of rj. Lines repre- 
senting the M—cr relation for constant values of rj are superimposed. 
The solid and dashed lines for which the power and intercepts of 
the M-a relation log 10 (M BH / M ) = B + <5 log 10 (cr/200 km s -1 ) 
are (S,B) = (4.0,8.1), (5.0,8.1), respectively. See the electronic 
edition of the Journal for a color version of this figure. 

the literature. Lines of constant rj for M-a relations 
log 10 (M B H/MQ) = J B + <51og 10 (cr/200kms~ 1 ), for which 
the power and intercepts are (5, B) = (4.0, 8.1), (5.0, 8.1), 
are superimposed in solid and dashed, respectively. Sim- 
ulation points are shown in blue and green with different 
shades of fill color denoting the different filling factors: 
fv = 0.42, 0.13, 0.052, 0.027 correspond to filled, dark 
grey, light grey, and white. Square, round, and triangular 
markers represent simulations for which fc m i n = 40 kpc -1 , 
20kpc~\ lOkpc -1 , respectively. 

The interpretation of the location of a point for the 
outflow of an observed galaxy on this plane is the same 
as that for the simulation points. If the point lies above 
a given line of constant 77, feedback is effective in that 
galaxy if that value 77 applied to that galaxy. For some 
galaxies we have calculated the value of rj = r/ \y S , so 
a direct comparison between that and the critical value 
?7 cr it marked by its location with respect to the lines of 
constant 77 can be made. We find that most galaxies with 
outflows listed in Table [3] lie below the line ?7 cr it > 10~ 3 . 
For most galaxies ?7obs > ?7crit, meaning that feedback 
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Figure 12. The absolute frequency N of the ratio of outflow 
velocity to the velocity dispersion, v ou t / cr^as predicted by the M— 
a relation for all radio galaxies in Table [3] for which the outflow 
velocity and SMBH mass are available. The 14 bins in v ou t/a each 
have a width of unity. The first orange bin counts the number 
of sources for which feedback is not effective, while the green bins 
count the sources in which feedback is effective. The darker colored 
bars denote outflow velocities of ionized gas and the lighter colored 
bars those measured for the neutral gas component in the galaxies. 
The unfilled bars show the combined distribution for all outflow 
velocity measurements; those above the darker and lighter bars, 
are, respectively, for the distributions for which the measurement of 
the neutral gas was preferentially take over that for the ionized gas, 
and vice versa. Overlaid are the values of r\ (right ordinate axis) 
for each binned galaxy. See the electronic edition of the Journal 
for a color version of this figure. 



associated with the outflows is effective in these galaxies. 

The effectiveness of feedback is better seen in terms of 
the distribution of the ratio of observed outflow velocity 
to velocity dispersion predicted by the M-a relation for 
the galaxies, v out /a. Here we assume an M-a relation 
for which the power and intercept are 5.3 and 8.2, respec- 
tively, values that give the best fit to elliptical galaxies, 
according to |Graham et al. | ([ 201 1 ^ . The distribution is 
shown as a histogram Fig. |12} Each of the 14 bins has 
a width of unity. Thus, the bin in orange to the left 
of v out /a = 1 counts the number of sources for which 
feedback is not effective (according to the criterion used 
in this paper), whereas the green bins to the right of 
f out /a = 1 count the sources for which feedback is effec- 
tive. Again we differentiate between observations of the 
ionized and neutral phases; the left, darker bars in a bin 
count the number for measurements of the outflow veloc- 
ities in ionized phase, and the right, lighter bars count 
those for the neutral phase. The unfilled bars trace the 
combined distribution of outflow velocities. The unfilled 
bars above the bars for the ionized (neutral) gas outflow 
velocity distribution show the combined distribution for 
which the measurement of the ionized (neutral) outflow 
velocity was preferentially taken over the measurement 
of the neutral (ionized) outflow velocity. Scattered on 
the same plot against the right ordinate axis are values 
for X] for eac h o f the binned items. 

From Fig.[l2]we find that the outflows in most galaxies 
have velocities that are a factor of a few greater than a 
velocity dispersion predicted by the M-a relation. Only 
in a few galaxies the outflows are weak in relation to the 
mass of the galaxy. There is no correlation between the 
values of r\ for galaxies in which feedback is effective and 
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in those in which feedback is not effective. In both do- 
mains, the scatter is large. However, most galaxies have 
large values of rj, with n > 10~ 2 . The locations of the 
points for the observed galaxies in Fig. |ll| in relation to 
those of our simulations suggest that the ISM distribu- 
tion in radio galaxies where fast outflows are observed is 
clumpy on quite large scales (i? c ,max > 25 pc) and possi- 
bly quite dense with a mean hot phase density of 1 cm~ 3 
within the inner 1 kpc. 

Because of the small sample size and large uncertain- 
ties associated with observationally measured jet pow- 
ers, outflow velocities, and SMBH masses, none of the 
conclusions above are definitive. The comparison be- 
tween observations and simulations presented here are 
merely a first step at understanding how effective feed- 
back is in radio galaxies at different redshift. A benefit 
of this method is that one may draw conclusions about 
the gas distribution of the hosts' ISM without direct spa- 
tially resolved observations. Conversely, if spatially re- 
solved observations of the hosts' ISM become possible 
with, e.g. the Atacama Large Millimeter/subm illimeter 
Array (ALMA, |Planesas||2011| jLons dalc 201 1, or the 
Square Kilometer Array (SKA, |Uodtrey et al.|2012| ), the 
parameter space of relevant simulations and the derived 
feedback efficiencies can be strongly constrained. A simi- 
lar comparative approach as demonstrated here between 
theory and observations could be made for radio quiet 
AGN that feature other modes of feedback. 

6.2. Other negative and positive feedback criteria 

The criterion for effective negative feedback adopted in 
this work that the warm-phase outflow velocities driven 
by AGN jets exceed the velocity dispersion predicted by 
the M-a relation for a galaxy is not unique. The cho- 
sen criterion, however, has the advantage that it is di- 
rectly relev ant to the evolution of the ISM in the bulge 
of galaxies ( Kormendy et al.||2009 ) , and that simulation 
and observation may be superimposed on the same plane 
spanned by velocity and jet power. One alternative is 
to use the condition that the material in the jet-driven 
outflow reaches the escape velocity and will be unbound 
from the potential of a galaxy. The fate of the outflowing 
gas is an interesting question in itself, and will need to be 
determined by simulations on scales of 10 - 1000 kpc sim- 
ulations. Furthermore, the fulfilment of either condition 
does not necessarily imply inhibited star-formation. 

Positive feedback is more difficult to quantify because 
the conditions for star formation are governed by a range 
of competing physics influencing the gravitational stabil- 
ity of clouds and the formation of dense cold cores within. 
They include external pressurization, hydrodynamically 
and conductively driven ablation, shocks driven into 
clouds, X-ray ionization, and molecular chemistry, and 
non-ideal MHD effects. Jet-induced star- format ion is 
thought to occur in gas-rich galaxies and pr oto-clusters 
at high redshifts ([Miley & De Breuck|[2008| , but there 
are als o examples in the nearby universe, e.g. in C entau- 
rus A jMould et al.|2000| , and Minkowski's object ([C roft 
eT^[l)05| , and 3C 285 ( |van Breugel fc DeyjllQ^ Tm' 
which the radio jet is implicated in shock- and pressure- 
trig gered star formatio n. 



sure bubble driven by the jet compresses the gas in the 
disc. We also observe an increase in the probability dis- 
tribution of dense gas in our simulations. Both external 
compression and (under-resolved) radiative shocks con- 
tribute, and to properly assess whether the a cloud would 
collapse to form stars, it is necessary to perform simu- 
lations with self-gravity. WB11 compared the mass of 
the clouds to the critical Bonnor-Ebert mass before and 
after the clouds were engulfed by the high-pressure jet 
plasma and concluded that the external pressurization 
places the clouds in the unstable regime, but also that 
the cloud ablation time-scales were short compared to 
the collapse timescales. A simil a r conc lusion was reached 
by Antonuccio-Delogu & Silk (2008) with simulations 
of a jet passin g near an isolated cloud. The fact that 



Gaibler et al. (2011) showed that the star-formation 
rate in a disc galaxy is enhanced a factor of 2 - 3 as a 
result of jet-ISM interactions, primarily because the pres- 



Gaibler et al.| (| 2011[ ) observe a marked increase in the 
star formation rate in the compact disk of molecular gas 
with high filling factor is consistent with the picture that 
cloud ablation is less impor tant the larg er the cloud com- 
ple xes are, as described i n |5.2|and { 5.3 The simulations 
by Sutherland fc Bicknell| (|2007) also demonstrate that 
dense gas distributed in a disc-like geometry couples less 
readily to the jet in the form of an outflow. It is clear 
that the consequences of AGN jet feedback (and other 
modes of feedback) depend as much on the multi-phase 
ISM properties of the galaxy during feedback as on the 
power of the central source. 

6.3. Cloud ablation 

The ablation and destruction of clouds is a difficult 
process to capture accurately in hydrodynamical simu- 
lations. There have been many studies investigating the 
destruction of clouds overrun by shocks (e.g. in super- 
nova remnants), or clouds embedded in a flow (e.g. a stel- 
lar wind) with numerical simulations, taking i nto account 
a variety of physical effects including coolin g ( Vietri et al. 
1997) iCooper et al.|2009] jYirak et al.|2010[ ) , thermal con- 
duction ([Orlando et al. 2005, 2008), structural inhomo- 
geneity (jXu'F'S tonc 1995; Poludncn ko et al |2002l|Naka- 
mura et al. 2006; Cooper et al. 2009), and magnetic fields 
( |Gregori (M mTTJI !)l)9[ [Shin et al ||20D8| ). Rayleigh- Taylor 
and Keivin-Helmholtz instabilities create a turbulent in- 
terface between the cloud surface and ambient flow where 
the mixing between the two phases occurs, but resolution 
limitations and the artificial viscosity due to numerical 
diffusion in hydrodynamic simulations clamp the spatial 
scales and energy scales th at can be captured, thereby 
diluting the mixing process ( TPittard et aLl[2009| . 

Since the comprehensive a nalytic and numerical work 
by Klein et al. (1990 1994), the fiducial minimum res- 
olution to capture the complete destruction of an adia- 
batic spherical uniform cloud in hydrodynamical simula- 
tions has been accepted to be around 100 cells. Radia- 
tive cooling, however, ra dically changes the destruction 
mech anism of the clouds ( Vietri et al.|1997 Cooper et al. 
2009 1 . In most astrophysical flows including those in our 
simulations, the cooling timescale is much shorter than 
the flow crossing time, and radiative shocks driven into 
clouds develop a thin, protective wall near the edge of the 
cloud boundary, which inhibits the Keivin-Helmholtz in- 
stability from rapidly destroying the cloud. | Cooper et al. 
(2009) showed that the cloud breaks up into long-lived 
filamentary cloudlets advected far downstream. 

While individual clouds in our simulations are under- 
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resolved, it is not clear by how much we are systemat- 
ically overestimating or underestimating cloud ablation. 
The fact th at we obtain a value of the ablation coefficient 
defined by ( Hartquist et al.p986 ), a, in Eqn. (261 of or- 
der 10, rather than I, from the predicted mean bubble 
density at late times of our simulations may indicate that 
we are overestimating the cloud ablation rate in our sim- 
ulations, perhaps because the stabilizing cooling layers 
behind the clouds are insufficiently resolved. 

The development of turbulence and the degree of frag- 
mentation during shock-cloud and shock-wind interac 



tions depend on the physical structure of the clouds ( Xu 
fc Stone|[T995l IPoludnenko et al.||2002l |Nakamura etal 



Cooper et al.||2009[) . |Nakamura et al.| |2006p r~for 
>le, find that cloud destruction is prolonged for 



2006 

example, find that cloud destruction is prolonged for 
smoother interfaces between c louds and the embed ded 
medium. On the other hand, Cooper et al.| ( |2009[ ) ob- 
serve that fractal clouds fragment taster than spheri- 
cal clouds. The large value of ablation coefficient may, 
therefore, also be a realistic consequence of the inhomo- 
geneous, fractal outlines of our clouds, which seed the 
Kelvin-Helmholtz instabilities and initially enhance the 
ablation rate. The ablation and entrainment rate seen in 
our simulations may also be higher because the flow, in 
which the clouds are embedded, is already rather turbu- 
lent as the jet streams percolate through the inter-cloud 
channels. Pressure variations and "buffeting" , particu- 
larly in the wake of the cloud may be efficient at extract- 
ing material from the back of the clouds and entr aining it 
into the tail streams. This effect was observed by [Pittard] 
|et al.| |2009[) to l e ad to faster cloud destruction. 

[ Pittard et al. (2010) directly compare the lifetimes 
of clouds in their simulations with those predicted by 
Eqn. (26) and find that, for clouds with low density con- 



trasts embedded in low Mach number flows the expres- 
sion underestimates the ablation rate by a factor of ~ 4, 
while for clouds with high density contrasts embedded in 
high Mach number flows the predicted values are a fac- 
tor of ~ 2 - 5 larger than found in the simulations. This 
may also partially account for the large value of a found 
in this work. 

The statistical distribution of clouds in our simula- 
tions fo r which fc m j n = 20 kpc - 1 , is identical to that 
used by Sutherland & Bicknell (2007), who performed 
a three level resolution study showing that the fractal 
features in their simulatio ns were at least partly cap- 
tured. As demons trated by Stone & Norman ( 1992[) and 



Klein et al.| ( 1994 ) for spherical adiabatic clouds7 Cooper 
et ai. "( 2009 ) found in their study of individual radiative 



fractal clouds embedded in a supersonic flow that the 
rate of fragmentation depended on the resolution of the 
cloud. However, no convergence was found at the high- 
est resolutions. A re cent resolution convergence study 



by Yirak et al. (2010), probing a broader range of reso- 
lutions of radiative shock-cloud interactions, found that, 
unlike adiabatic cases, the flow structure does not show 
signs of convergence at 100 cells per clump; convergence 
may only formally be reached when the cooling length is 
well resolved. This is highly impractical in global simula- 
tions such as those conduc ted in this work. O ne possible 
improvement discussed by Yirak et al. ( |2010 |, although 
complicated, involves a careful use of adaptive mesh re- 
finement of the cooling layers in the radiative shocks. 
Another possibility is to implement a subgrid treatment 



of turbulence that leads to better con vergence with in- 
creased resolution ( Pittard et al.||200"9 ). 

Given the difficulty m ensuring sufficient resolution to 
reach flow convergence and the complexity of including 
all the physical effects that may modify the cloud de- 
struction timescale in opposing ways, the setup of our 
simulations, despite limited resolution across one cloud 
are justified as a first step to model AGN jet feedback in 
fully three-dimensional hydrodynamic simulations. 

6.4. Long-term evolution 

The domain extent in our simulations covered the cen- 
tral 1 kpc 3 of a gas rich radio galaxy, within which the 
coupling of radio jet to the dense gas is strongest. The 
fate of the outflowing gas on scales larger than 10 kpc 
will influence the accretion and star-formation rates of 
the galaxy at later times and is also relevant to theo- 
ries of the enrichment of the intracluster medium (ICM). 
An example in which an AGN radio jet may be directly 
responsible for carrying enrich ed material into the ICM 
was id entified in 4C+44.16 by |Hlavacek-Larrondo et al.| 
(20111. But the fraction of gas that is unbound from 
the galaxy cannot be predicted from our simulations and 
separate simulations on larger scales with the inclusion of 
a gravitational field need to be performed. In this prob- 
lem, non-uniformity of the dense gas distribution and 
asymmetries in the energy injection influence the require- 
ments to unbind gas from a d eep gravitational potential 
(Bland-Hawthorn et al. 2011). An inhomogeneous ISM 
and off-center energy injection reduce the fraction of gas 
that can reach escape velocities. 

In a gravitating environment, buoyantly rising jet- 
inflat ed bubbles and associa ted buoyancy-driven instabil- 
ities ( Balbus & Soker 1989[ ) will influence the mi xing rate 
and evolution of the ISM ( jBriiggen et al.||2002[ ). Large- 
scale simulations will also aid comparisons with observa- 
tions of the radio jet and the outflowing gas, which have 
limited resolution on kpc scales. 

The balance between heating and cooling in the I CM of 
cool-core clusters can also be tested. For example, |Gas-| 
pari et al. (2012) performed grid-based hydrodynamic 
simulations on these scale and found that a feedback cy- 
cle involving a heavy, slow jet reproduced the typical 
observed entropy profiles around cool core cluster cen- 
tral galaxies. Such heavy, slow jets may be the result 
of entrainment and mass-loading that jets experience in 
the core of galaxies though interactions similar to those 
shown in our simulations. Gaspari et al. also demon- 
strated the pervasiveness of a two-phase ISM owing to 
cold phase gas condensing out of the hot phase through 
the thermal instability and fuelling the AGN. 

7. CONCLUSIONS 

We have conducted a total of 29 3D grid-based hydro- 
dynamical simulations of AGN jets interacting with a 
fractal two-phase ISM. The simulations cover jet powers 
of Pjet = 10 43 ergs -1 - 10 46 ergs -1 and a range of differ- 
ent ISM parameters characterized by the density of the 
hot phase, the filling factor of the warm phase (clouds), 
and the maximum cloud sizes in the fractal distribution 
of the warm phase. The simulations are applicable to the 
early and intermediate phases of radio-mode feedback at 
high redshifts, which often involve massive, evolving gas- 
rich (proto)-galaxies. The sufficiently broadly sampled 
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parameter space allowed us to study the dependence of 
the negative feedback efficiency, as measured by the max- 
imum outflow velocities reached through the jet-ISM in- 
teractions, on the filling factor and the maximum clouds 
sizes. We also identified the precise physics of the mo- 
mentum and energy coupling that leads to rapid acceler- 
ation of the warm phase. Finally, we undertook a com- 
parison between recent observations of outflows in radio 
galaxies and our simulations on the plane defined by out- 
flow velocity and jet power. 

Fifteen new simulations were conducted to supplement 
the simulation data by WBlf with results for feedback 
efficiencies for filling factors of 0.052 and 0.027, and 
those with maximum cloud sizes of -R c ,max = 10 pc and 
-Rc.max = 50 pc. The lower filling factor runs and runs 
for which the maximum cloud sizes were i? c ,max = 10 pc 
were conducted at a resolution of 1 pc per cell, twice the 
resolution of the simulations by WB11 within the same 
box size. The principal findings from these new simula- 
tions are the following: 

f. The main conclusions reached in WB11 remain 
unchanged: Feedback is effective in systems in 
which the jet power is in the range Pj 0t = 10 43 
- 10 46 ergs -1 and r\ > ?7crit, where r\ is the Edding- 
ton ratio of the jet, 77 = Pj e t/L e dd- The critical 
Eddington ratio, r) CI it, is determined by the crite- 
rion that the velocity of the outflow driven by a 
jet of power Pj et exceed the velocity dispersion ex- 
pected from the M-a relation of a galaxy for which 
the Eddington luminosity is L e dd- For reasonable 
values of the ISM parameters, feedback ceases to 
be efficient for Eddington ratios of r\ < 10~ 4 . Due 
to a large sustained mechanical advantage in these 
systems, the fraction of jet energy transferred to 
the warm phase is ~ 0.1 - 0.4. 

2. The dependence of the feedback efficiency on fill- 
ing factor found by WB11 reverses for filling factors 
fv 0.1 in the sense that a smaller filling factor 
leads to higher feedback efficiencies. In general the 
dependence on filling factor is weak, however. The 
reversal occurs because of a shift in the balance 
of two opposing effects as the filling factor is de- 
creased: the increase in surface area of the clouds 
exposed to ablation relative to their volume, and 
the increase in the volume between clouds through 
which the jet plasma may flood. The latter de- 
creases the feedback efficiency while the former in- 
creases the feedback efficiency and dominates below 
fv ~ 0.1. 

3. For a given filling factor, the feedback efficiency 
is higher the smaller the maximum cloud sizes 
in the fractal distribution. Galaxies containing 
cloud complexes initially larger than 50 pc (/c m j n = 
10 kpc~ ) require jets with Eddington ratios rj > 
10~ 2 for efficient negative feedback. Pressure trig- 
gered star-formation may be expected for large 
clouds (> 50 pc), while clouds smaller than ~ 10 pc 
(fcmin = 40kpc -1 ) are unlikely to survive signif- 
icant ablation. The dependence of feedback effi- 
ciency on cloud size is much stronger than that on 
filling factor and scales nearly linearly with fc m i n . 



This is the result of the linear relation between the 
size of clouds and the surface area of the clouds ex- 
posed to ablation relative to their volume. We in- 
troduce the concept of an interaction depth for jet- 
cloud interactions analogous to optical depth. The 
interaction depth increases linearly with smaller 
cloud sizes leading to a linear increase in outflow 
velocity with fc m i n . 

4. A comparison between the global dynamics of the 
outflowing warm phase material with that of an 
energy-driven bubble shows that outflows approach 
the energy-driven limit in cases for which the fill- 
ing factor is low or the maximum sizes of clouds 
is large. For a given filling factor, the dispersal of 
clouds is higher if clouds are smaller. Conversely, a 
jet-driven bubble impacting and engulfing a dis- 
tribution of large clouds may lead to pressure- 
triggered star-formation. Thus, the size distribu- 
tion of clouds strongly influences whether feedback 
is negative or positive. Considering the relatively 
weak dependence of bubble expansion rate on cloud 
sizes, we argue that, if conditions in the ISM of ra- 
dio galaxies in cosmological simulations are in the 
regime of fv < 0.1 and -R c ,max ^5 25 pc, an energy- 
driven sub-grid implementation of (negative) AGN 
radio mode feedback is justified. In this regime, 
feedback in a single phase medium is a good ap- 
proximation to feedback in a two-phase medium 
because the warm phase material embedded in the 
energy-driven bubble is accelerated to speeds com- 
parable to the bubble expansion speed within the 
dynamical time of the bubble. 

5. We find that a simple theory of clouds embedded 
in a fully thermalized energy-driven bubble does 
not provide sufficient ram pressure to accelerate 
the clouds to velocities observed in high redshift 
radio galaxies. Instead, the momentum transfer is 
provided by the ram-pressure of the partially ther- 
malized streams of jet plasma flooding through the 
inter-cloud channels, which have turbulcntly en- 
trained ambient hot gas (initially external to the 
expanding bubble) and are mass-loaded with ab- 
lated warm cloud material. Initially the chan- 
nel flows carry turbulently entrained shocked hot- 
phase material at velocities of ~ 10 kms , and 
at later times, the channel flows are dominated by 
material ablated from clouds. The acceleration ef- 
ficiency is highest in the former stage, while mass- 
loading from clouds reduces the acceleration effi- 
ciency somewhat. This mechanism, which is rem- 
iniscent of the two-stage feedbac k mechanism pro- 
posed by Hopkins & Elvis ( |2010[ ), is capable of ac- 
celerating the clouds to velocities of 100s to 1000s 
kms -1 within the dynamical time of the bubble. 

6. The observed outflows of neutral and ionized ma- 
terial in most of the radio galaxies in our sam- 
ple compiled from the literature are fast enough to 
cause substantial velocity dispersions in the host. 
The critical jet Eddington ratios in most sources is 
?7crit ^ 10 -3 , but the observed jet Eddington ratios 
are predominantly rj > 10 -2 . By comparing the jet 
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powers and outflow speeds obtained from the ob- 
servations with those found in our simulations, we 
tentatively infer that the ISM of the radio galaxies 
in which outflows are observed is clumpy on quite 
large scales (i? c ,max > 25 pc) and possibly quite 
dense with a mean (hot phase) density of 1 cm~ 3 
within the inner 1 kpc. 

7. We explain some of the radio-mor phological char 



acteris tics of 3C 48 discussed by |Stockton et al. 
(20071 with a synthetic radio image of one of our 
simulations, including the jet collimation within 
the extended O in emission region, and the deflec- 
tion and expansion of the jet lobe beyond the emis- 
sion region. 
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